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Abstract 

Kolmogorov-Arnold-Moser  (KAM)  Theory  states  that  a  lightly  perturbed, 
conservative,  dynamical  system  will  exhibit  lasting  quasi-periodic  motion  on  an  invariant 
torus.  Its  application  to  purely  deterministic  orbits  has  revealed  exquisite  accuracy  limited 
only  by  machine  precision.  The  theory  is  extended  with  new  mathematical  techniques  for 
determining  and  predicting  stochastic  orbits  for  Earth  satellite  systems.  The  linearized 
equations  of  motion  are  developed  and  a  least  squares  estimating  environment  is 
pioneered  to  fit  observation  data  from  the  International  Space  Station  to  a  phase  space 
trajectory  that  exhibits  drifting  toroidal  motion  over  a  dense  continuum  of  adjacent  tori. 
The  dynamics  near  the  reference  torus  can  be  modeled  with  time-varying  torus 
parameters  that  preserve  both  deterministic  and  stochastic  effects.  These  parameters  were 
shown  to  predict  orbits  for  days  into  the  future  without  tracking  updates — a  vast 
improvement  over  classical  methods  of  orbit  propagation  that  require  routine  updates. 
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STOCHASTIC  ORBIT  PREDICTION  USING  RAM  TORI 


I.  Introduction 


1.1  A  Silent  Crisis 

Students  from  the  Sewickley  Academy  packed  tightly  into  a  small  room  at  the 
Carnegie  Science  Center  in  Pittsburg  on  the  morning  of  March  12,  2009,  for  an  amateur 
radio  downlink  from  the  International  Space  Station  (ISS).  One  of  the  students  waited 
with  nervous  excitement  to  ask  about  “new  objects  being  sent  into  space  to  be  attached  to 
the  ISS,”  but  ironically  he  was  about  to  get  a  rapid  lesson  about  objects  sent  into  space 
that  aren’t  attached  to  anything  useful  at  all,  i.e.  frangible  nut  fragments,  drifting  tool 
bags,  defunct  satellites,  loitering  rocket  stages  the  size  of  school  buses  -  junk  in  general 
[1;  2]. 

Just  moments  before  the  scheduled  downlink,  Mission  Control  Center  (MCC) 
Houston’s  Capsule  Communicator  (CAPCOM)  Kathy  Bolt  delivered  a  chilling  message 
to  Expedition  18  Commander  Michael  Fincke: 

We  have  information  about  a  RED  conjunction.  The  information  came  in  late  and 
with  the  uncertainty  of  it,  we  are  wanting  to  take  a  conservative  approach.  It ’s  a  low 
probability  of  a  hit;  however,  the  object  is  rather  large  based  on  what  we  can  track 
and  if  it  does  happen  to  hit  the  ISS,  we  ’re  talking  about  only  a  10  minute  reserve  time. 
[3] 
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The  13-cm  diameter  object,  cataloged  as  “25090  PAM-D,”  was  a  tumbling  yo  weight  that 
unwound  from  a  spent  payload  assist  module  (PAM-D)  during  the  launch  of  Navstar  2A- 


11  in  1993  [4;  5], 


As  a  precaution  against 
depressurization  from  a  collision, 
the  three  person  crew  took  refuge  in 
the  Russian  Soyuz  “lifeboat”  shown 
in  Figure  1.  Lucky  for  them,  the 
object  passed  without  incident,  but 
this  was  a  close  call  that  could  have 


been  avoided.  Normally  the  ISS  is 

Figure  1.  Russian  Soyuz  docked  to  the  ISS. 
maneuvered  out  of  harm’s  way.  In  Credit:  NASA 

the  eight  years  prior  to  this  event,  the  ISS  maneuvered  eight  times  to  avoid  space  debris 

[6].  No  maneuver  was  coordinated  this  time  because  25090  PAM-D  was  classified  as  a 

low  concern  42  hours  prior  to  the  time  of  closest  approach.  NASA  claimed  to  have  “good 

tracking  and  a  miss  outside  the  notification  box”  [4],  but  after  the  maneuver  window 

passed,  it  was  discovered  that  the  solar  radiation  pressure  in  the  tracking  model  was 

wrong  due  to  the  low  perigee  of  the  object.  This  resulted  in  a  smaller  radial  miss  distance. 

NASA  scrambled  to  obtain  a  higher  certainty  of  the  object’s  path,  but  the  highly  elliptical 

and  low  perigee  orbit  made  it  difficult  to  get  a  good  track.  The  best  models  showed  a 

stable  and  well-behaved  covariance,  but  it  was  much  larger  than  typically  observed  near 

the  time  of  conjunction  [4], 
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This  wasn’t  the  first  time  NASA  encountered  such  shortcomings,  nor  would  it  be  the 
last.  Just  eleven  days  later,  the  ISS  had  to  be  maneuvered  to  avoid  object  26264,  part  of  a 
Chinese  CZ-4  rocket  launched  in  1999  that  broke  up  in  March  2000  [7].  These  close  calls 
were  warning  shots  fired  across  the  bow  of  the  ISS.  It  made  painstakingly  clear  the 
limitations  of  our  tracking  resources,  capabilities,  and  the  certainty  of  orbit  predictions, 
all  of  which  will  habitually  haunt  us  in  the  absence  of  immediate  and  pervasive 
innovations  to  the  way  we  conduct  business. 

1.2  Purpose 

The  subject  of  this  thesis  will  explore  a  new  method  for  predicting  orbits  that  may 
improve  our  space  situational  awareness  and  allow  us  to  track  more  objects  with  greater 
accuracy  than  we  can  today.  Doing  so  would  potentially  enhance  space  operations  to 
more  easily  avoid  debris  headed  our  way.  Unfortunately,  even  the  most  accurate  method 
of  orbit  determination  cannot  eliminate  the  excessive  catalogue  of  junk  orbiting  the  earth. 
That  makes  it  all  the  more  essential  that  we  have  the  best  orbit  predictions  possible  to 
forecast  collisions  and  take  whatever  actions  are  available  to  reduce  the  impact  on  other 
operations. 

1.3  Implications 

Since  the  dawn  of  the  space  age,  the  United  States  and  other  space  faring  nations  have 
lobbed  manmade  objects  into  orbit  under  the  premise  that  space  is  so  big,  it  is  highly 
unlikely  that  two  random  objects  will  collide  in  a  limited  time  span.  This  paradigm  - 
called  the  Big  Sky  Theory  -  had  left  us  more  concerned  about  junk  falling  out  of  the  sky. 
Routine  encounters  with  space  junk  reentering  the  atmosphere  supported  this  early 


3 


outlook.  Take  for  example  the  Lan  Chile  Airbus  A340  which  was  travelling  between 
Santiago,  Chile,  and  Auckland,  New  Zealand  on  March  27,  2007  [8].  The  pilot  was 
enjoying  the  calm  night  air  over  the  South  Pacific  Ocean  when,  out  of  nowhere,  flaming 
balls  of  space  junk  went  hurtling  past  the  plane.  The  debris  was  later  identified  as  a 
reentering  Russian  Progress  23P  cargo  freighter  for  the  ISS.  Similar  incidents  had  a  more 
direct  “impact”  on  our  psyche.  Skylab’s  premature  reentry  in  1979  wiped  out  a  cow  in  the 
Australian  Outback  and  in  1997  an  Oklahoma  woman  was  struck  in  the  shoulder  by 
charred  material  from  a  Delta  II  rocket  that  the  U.S.  Air  Force  launched  the  previous 
year. 

With  the  threat  of  space  debris  coming  from  all  directions,  U.S.  Space  Command’s 
Space  Surveillance  Network  (SSN)  today  tracks  over  19,000  objects  in  earth  orbit 
roughly  the  size  of  a  baseball  or  larger.  Even  as  they  prioritize  tracking  assets  to  protect 
human  spaceflight  missions  and  high  priority  military  and  civilian  satellites,  prediction 
accuracy  is  limited  and  hundreds  of  thousands  of  smaller  artificial  debris,  from  paint 
flecks  to  solid  rocket  fuel  slag,  are  untracked.  The  smaller  objects  can  pose  just  as  serious 
of  a  threat  in  the  space  environment,  but  their  sheer  number  and  size  make  tracking  these 
objects  very  difficult.  With  a  lower  mass-to-area  ratio,  air  drag  is  a  more  dominant  force, 
tending  to  degrade  the  orbit  much  too  rapidly  for  sustainable  predictions. 

As  the  population  density  of  satellites  and  debris  increases,  the  finite  probability  of 
two  objects  colliding  also  increases.  When  probability  becomes  reality,  the  density 
skyrockets  further  with  the  production  of  new  fragments  that  pose  new  threats  to  existing 
satellites.  For  example,  in  February  2009  when  the  Russian  Cosmos  2251  collided  with  a 


4 


U.S.  Iridium  satellite  500  miles  above  Siberia,  it  resulted  in  a  vast  cloud  of  debris  adding 
to  the  already  crowded  low-altitude  orbital  catalogue.  More  than  2,000  objects  being 
tracked  today  are  a  result  of  the  Cosmos-Iridium  collision  alone. 

In  1978,  NASA  scientists  Donald  Kessler  and  Burton  Cour-Palais  predicted  that  a 
continued  population  growth  of  satellites  coupled  with  seed  collisions,  like  the  one 
described  above,  would  fuel  a  cascade  of  collisions.  “The  result  would  be  an  exponential 
increase  in  the  number  of  objects  with  time,  creating  a  belt  of  debris  around  the  earth 
[9].”  This  is  very  much  like  the  natural  planet  forming  process  that  relies  on  a  cascade  of 
collisions  beginning  with  larger  objects  and,  over  time,  shifting  to  smaller  objects  that  are 
greater  in  number.  With  every  collision,  the  smashed  remnants  lose  their  inclination  and 
eccentricity  until  eventually  a  cloud  of  dust  orbits  in  a  ring  about  the  equator  of  the 
central  body  (like  the  rings  of  Saturn),  or  coalesce  into  a  moon  or  planet.  Hannes  Alfven 
famously  described  this  coalescing  of  debris  using  apples  in  a  spacecraft  [10]. 

Later  studies  done  by  Kessler  identified  critical  population  densities  for  unstable  and 
runaway  environments  in  which  collisional  cascading  occurs  [11;  12],  He  surmised  that  a 
threshold  value  for  instability  exists  in  an  environment  in  which  the  fragment  population 
will  inevitably  increase  as  a  result  of  random  collisions  dominated  by  overpopulated 
“intact”  objects.  Given  enough  time  and  no  additions  to  the  intact  population  (only 
subtractions  by  way  of  collision),  equilibrium  will  be  achieved  in  which  the  number  of 
fragments  generated  by  random  collisions  matches  the  number  of  fragments  eliminated 
by  natural  processes  (like  air  drag).  Increasing  the  intact  orbital  population  beyond  its 
final  state  at  equilibrium  would  upset  the  system  and  the  production  rate  of  debris  would 
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exceed  the  decay  rate  of  debris.  Under  the  scenario  that  the  intact  population  is 
continually  increased  or  held  steady  above  the  equilibrium  threshold,  the  number  of 
fragments  would  increase  toward  infinity.  Kessler  describes  this  as  a  “runaway 
environment”  [13]. 

The  NASA  Orbital  Debris  Program  Office’s  LEGEND  model  for  studying  the  orbital 
debris  environment  between  200  and  50,000  km  altitude  currently  shows  regions  between 
600  and  1,700  km  altitude  that  are  already  beyond  the  critical  density  [14].  Figure  2 
shows  a  200  year  forecast  for  the  altitude  band  between  900  and  1,000  km.  With  no  new 
launches  contributing  to  the  current  intact  population  of  roughly  500  objects  in  that  band, 
the  prediction  shows  a  runaway  environment.  Similar  conditions  are  found  near  1,400  km 

[13]. 


Year 


Figure  2.  LEGEND  prediction  shows  runaway  environment  between  900  and  1,000  km 
assuming  no  new  additions  to  the  intact  population  after  2004. 
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Earlier  last  year,  Kessler  warned  that  if  the  aerospace  community  continues  to 
conduct  business  as  usual — without  international  adherence  to  post  mission  disposal 
guidelines — and  if  we  refuse  to  innovate  with  active  debris  removal  programs, 
catastrophic  collisions  will  become  more  frequent.  This  could  wreak  havoc  on  the  global 
economy  by  disrupting  commerce,  communications,  agriculture,  transportation,  health 
services,  education,  energy,  and  the  environment,  all  of  which  have  become  dependent  on 
a  diverse  network  of  space-based  systems  that  detect  and  transmit  information.  Just  as 
importantly,  it  could  impair  the  United  States’  ability  to  defend  its  interests  at  home  and 
abroad. 

With  so  much  at  stake,  Congress  established  the  largest  study  ever  conducted  on  U.S. 
space  management  and  operations  in  the  National  Defense  Authorization  Act  for  Fiscal 
Year  2000.  The  Commission  to  Assess  United  States  National  Security  Space 
Management  and  Organization,  chaired  by  the  Honorable  Donald  H.  Rumsfeld,  released 
its  findings  on  January  11,  2001,  with  a  slew  of  warnings  about  space  situational 
awareness.  “To  use  space  effectively  and  to  protect  against  threats  that  may  originate 
from  it,  the  U.S.  must  be  able  to  identify  and  track  much  smaller  objects  in  space  than  it 
can  track  today”  [15:31],  The  report  specifically  recommended  technological 
improvements  to  the  Space  Surveillance  Network  by  means  of  electro-optical  and  radar 
systems,  but  beyond  hardware  solutions  that  unveil  what  is  happening  this  instant,  the 
U.S.  strategy  must  also  include  proactive  measures  to  avoid  future  threats.  That  is  to  say, 
observing  precisely  what  struck  the  ISS  or  a  critical  satellite  is  less  important  than 
forecasting  that  something  will  strike  the  ISS  or  a  critical  satellite. 
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As  we  saw  earlier  with  the  close  call  at  the  ISS,  uncertainty  about  the  debris  path  left 
NASA  scrambling  for  better  tracking  data  to  corroborate  the  accuracy  of  already  limited 
predictions.  If  we  can  predict  orbits  with  greater  accuracy,  active  collision  avoidance 
(ACA)  systems  can  serve  as  a  partial  solution  for  those  objects  that  have  the  means  to 
maneuver,  such  as  the  ISS.  But  even  for  those  that  do  not,  there  are  still  benefits  to  be  had 
by  targeting  objects  for  new  disposal  initiatives  that  seek  to  remove  the  most  dangerous 
debris  from  LEO.  If  nothing  else,  an  expanded  catalogue  of  improved  predictions 
provides  a  crucial  level  of  situational  awareness  and  gives  us  the  ability  to  plan  rather 
than  react  and  potentially  avert  a  silent  crisis. 

1.4  State-of-the-Art 

The  current  state-of-the-art  for  orbit  propagation  can  be  divided  into  three  main 
branches:  numerical  (also  called  special  perturbations),  analytical  (also  called  general 
perturbations),  and  semianalytical  techniques  (a  combination  of  the  other  two).  The  most 
accurate  method  is  the  numerical  technique — generally  known  as  Cowell’s  method — that 
performs  a  direct  numerical  integration  of  the  equations  of  motion  with  expanded,  high- 
order  disturbing  forces.  The  problem  with  this  method  is  that  it  requires  very  small  time 
steps  to  benefit  from  the  precision  of  high  order  models  and  is  therefore  computationally 
expensive.  Even  as  computer  resources  have  become  pervasive,  the  sheer  number  of 
objects  that  must  be  tracked  makes  this  method  difficult  at  this  time.  The  analytical 
techniques  which  include  Simplified  General  Perturbations  (SGP),  its  popular  variant 
SGP4,  and  Brouwer-Lyddane  use  truncated  force  models  and  averaging  techniques  to 
simplify  the  calculations.  The  benefit  of  analytical  techniques  is  their  speed,  but  they 
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suffer  from  immediate  errors  at  epoch  on  the  order  of  1  km  due  to  the  simplifications 
[16].  Semianalytical  techniques  such  as  Draper  Semianalytic  Satellite  Theory  (DSST)  are 
hybrid  variants  of  the  previous  techniques  for  relatively  fast  and  accurate  results. 
Unfortunately,  all  three  methods  leave  a  lot  to  be  desired  since  even  the  most  accurate 
Cowell’s  method  can  become  suspect  for  low  earth  orbits  within  a  24  hour  period.  This 
lack  of  accuracy  keeps  the  Space  Surveillance  Network  (SSN)  overly  tasked  detecting, 
tracking  and  cataloging  as  many  objects  as  possible  -  barely  20,000  objects  out  of 
hundreds  of  thousands. 

Here  lies  the  inspiration  for  this  thesis:  what  new,  innovative  orbital  prediction 
techniques  can  we  validate  and  apply  to  accurately  track  more  objects  faster,  cheaper  and 
better? 

1.5  Problem  Statement  and  Approach 

This  thesis  will  introduce  new  mathematical  techniques  for  determining  and 
predicting  stochastic  orbits  based  on  the  work  of  William  Wiesel  [17;  18;  19]  using 
Kolmogorov-Amold-Moser  (KAM)  Theory  in  a  least  squares  estimation  environment. 
KAM  theory  eliminates  special  and  general  perturbation  paradigms  by  representing  orbits 
on  tori  in  6-dimensional  phase  space  rather  than  the  standard  rectangular,  inertial 
coordinates.  Wiesel  previously  showed  that,  when  perturbed  only  by  conserved  forces 
from  a  20x20  EGM96  gravity  field,  KAM  tori  predict  orbits  to  within  meters  of 
numerically  integrated  orbits  for  up  to  a  decade  [20]. 

With  such  high  accuracy  as  the  baseline  for  deterministic  orbits,  the  only  remaining 
hurdle  for  KAM  theory’s  application  to  low  earth  orbits  are  stochastic  effects  (air  drag, 
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lunar/solar  point  mass  effects,  and  other  non-conservative  perturbations).  Typically  for 
low  earth  orbits,  the  state-of-the-art  numerical  methods  diverge  from  the  actual  orbit 
within  a  couple  days  and  must  be  recomputed  from  new  initial  conditions,  whereas  this 
author  claims  that  the  baseline  KAM  torus  can  be  used  to  predict  the  actual  orbit  for  a 
much  longer  period  so  long  as  it  contains  routine  phase  angle  updates  and  stochastic 
offsets  to  the  reference  torus  (performed  in  this  research  at  every  time-step).  The  duration 
and  level  of  accuracy  of  this  so-called  “motion  near  a  reference  torus”  is  the  topic  of  this 
research. 

Previous  work  by  Little  attempted  to  fit  actual  observation  data  from  NASA’s 
oceanography  satellite,  Jason- 1,  and  the  Gravity  Recovery  and  Climate  Experiment 
(GRACE)  satellite  to  reference  tori  that  he  constructed  from  purely  deterministic  orbits 
[21].  His  results  for  Jason-1  produced  residuals  constrained  to  less  than  1  km  for  30  days 
which  is  not  surprising.  KAM  theory  states  that  a  lightly  perturbed,  nearly-integrable 
Hamiltonian  system  will  lie  on  a  quasi-periodic  trajectory  in  the  phase  space  of  a  torus. 
Since  the  height  of  perigee  for  Jason- 1  is  near  1,330  km — where  air  drag  no  longer 
dominates  among  perturbing  agencies  in  low  earth  orbit  (LEO)  [22:271] — the  system  is 
only  lightly  perturbed  and  is  governed  mostly  according  to  KAM  theory.  However,  the 
results  from  GRACE  depicted  a  much  worse  correlation  between  the  observation  data 
and  the  reference  torus  with  residuals  as  large  as  90  km.  This  is  due  to  frequent  thruster 
firings  and  air  drag  which  is  more  prominent  than  that  experienced  by  Jason- 1  since 
GRACE’S  nearly  circular  polar  orbit  has  a  height  of  less  than  500  km.  The  larger 
perturbations  cause  the  satellite  to  wander  off  the  reference  torus,  but  since  the  residuals 
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were  limited,  it  is  logical  to  suppose  the  satellite  was  constrained  to  motion  near  the 
reference  torus. 

Contingent  that  a  reference  torus  does  exist  for  a  non-chaotic  orbit,  Wiesel  and  this 
author  are  investigating  the  possibility  that  a  dense  continuum  of  flexing,  persisting  tori 
will  also  exist  nearby,  such  that  the  satellite  may  drift  off  the  initial  reference  torus 
(generated  only  by  Earth’s  gravity  field)  and  circumnavigate  neighboring  tori  when 
stochastic  perturbations  exceed  a  threshold  value.  This  must  be  true  for  KAM  theory  to 
accurately  apply  to  all  LEO  orbits,  especially  those  where  air  drag  is  most  prevalent. 

1.6  Research  Focus 

To  explore  the  concept  of  motion  near  a  reference  torus,  the  following  questions  and 
issues  will  be  addressed: 

1)  What  are  the  linearized  equations  of  motion  near  the  reference  torus? 

2)  What  classical  estimation  techniques  can  be  applied  to  fit  observation  data  to  a 
continuum  of  tori  near  the  reference  torus? 

3)  Given  a  good  fit  from  the  previous  question,  can  the  torus  be  used  to  generate 
stochastic  predictions  and  for  how  long? 

1.7  Preview 

All  three  research  focuses  will  be  successfully  answered  in  this  document.  Using  the 
ISS  orbit  as  a  test  case,  it  will  be  shown  that  a  single  “flexing”  torus  can  provide 
stochastic  predictions  with  constrained  residuals  for  a  yet-to-be-determined  time  limit 
(currently  more  than  18  days).  The  low  eccentricity  of  the  ISS  orbit  will  induce  obstacles 
that  limit  the  accuracy  of  these  predictions  to  RMS  residuals  near  2  km.  It  will  also  be 
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shown  that  the  torus  is  very  sensitive  to  drag  changes  in  LEO  orbits  and,  in  fact,  can  be 
used  to  identify  weak  variations  in  the  thermosphere  density  due  to  geomagnetic  activity. 
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II.  Background 


The  technical  merits  of  KAM  theory  are  anchored  firmly  in  Hamiltonian  physics 
which  necessitates  a  review  of  relevant  analytical  mechanics  in  this  chapter.  However,  as 
LEO  orbits  are  explored  in  this  research,  Keplarian  motion  (anchored  squarely  in 
Newtonian  physics)  must  be  relied  upon  more  heavily  for  understanding  motion  near  a 
KAM  torus.  This  necessitates  a  modest  review  of  orbital  mechanics  up  front.  Of  course, 
most  fundamental  to  this  research  is  KAM  theory,  itself.  Thus,  we  conclude  this  chapter 
with  a  look  at  the  most  fundamental  work  of  Kolmogorov,  Arnold  and  Moser. 

2.1  Orbital  Mechanics 

This  paper  is  intended  to  be  comprehensive  enough  to  avoid  routine  cross-references, 
but  in  order  to  stand  alone,  a  significant  review  and  derivation  of  the  two-body  problem  is 
necessitated  here.  The  two-body  problem  is  required  to  linearize  the  dynamics  of  motion 
near  Earth  satellite  KAM  tori  as  will  be  discussed  in  §3.2.  In  addition  to  deriving  the 
classical  orbital  elements,  a  number  of  other  related  topics  are  discussed  in  this  section, 
including  coordinate  systems,  perturbations,  and  Lagrange’s  Planetary  Equations — all  of 
which  are  relevant  or  essential  to  the  science  of  orbit  determination  using  KAM  tori. 

2.1.1  Keplerian  Motion 

In  1687,  Sir.  Isaac  Newton  published  The  Mathematical  Principles  of  Natural 
Philosophy,  more  commonly  known  as  the  Principia,  which  not  only  introduced  his 
famous  three  laws  of  motion,  but  also  the  law  of  universal  gravitation.  He  stated  that  any 
two  bodies  in  a  system  attract  one  another  with  a  force  proportional  to  the  product  of  their 
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masses  and  inversely  proportional  to  the  square  of  the  distance  between  them  [23:4],  For 
a  system  of  n  particles  given  by  mi,  m2,  mn,  the  vector  force  acting  on  mi  by  each 
particle  is  expressed  as  [23:7]: 


Fi  = 


ZTTlj 
7  —  1  Jl 


Jl 


j  =  1 
j*i 


(i) 


where  G  is  the  universal  gravitational  constant  approximated  as  6.672  x 
10"11  Nm2/kg2. 


The  position  and  velocity  vectors  for  the  center  of  mass  of  each  particle  is  given  by: 


r,  = 


Xi 

Vi 

LZi 


(2) 


Vi  = 


dxi  /dt 
dyi/dt 
dzi/ dt 


(3) 


The  force  acts  along  the  vector  [23:7], 

rn  =  r,  -  rs  (4) 

and  its  magnitude  is  expressed  as: 


For  a  simplified  case  with  no  external  forces  such  as  drag,  thrust,  solar  radiation 
pressure,  or  other  perturbations,  the  gravitational  force  can  also  be  formulated  in  the 
inertial  frame  using  Newton’s  second  law  of  motion  [23:8]. 


—  d  _  dVi  _  drill 

Fi  =  -T7  (miVi)  =  mi  —  +vi  — - 
dt  dt  dt 


(6) 
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Since  we  assume  the  body  is  not  expelling  mass  to  generate  thrust,  the  last  term  in  (6) 
is  zero. 

Setting  (6)  equal  to  (1)  and  dividing  both  sides  by  m;  gives  the  general  equation  of 
motion  for  the  ith  body.  A  total  of  n  second-order  vector  differential  equations  can  be 
formed  according  to  [23:9]: 

dt2  2-i  r.f  ^  (7) 

7=1  Jl 
j*i 


In  a  simplified  two-body  system,  the  two  resulting  second  order  differential  equations 
are  [23:9]: 


d2r1 

dt2 


m?  _  _ 

G  —  (r2  -  rx) 


'21 


d2r2 

~dtr 


mi  _ 

G  —  (rx  -r2) 


f  12 


(8a) 

(8b) 


From  (4)  we  know  that  [23:10]: 

d2r  12  d2r2  d2ri 
dt2  =~dir~~dtr 


Substituting  (8a)  and  (8b)  into  (9),  produces  the  simplified  Keplerian  (two-body) 
equation  of  motion  that  describes  the  acceleration  of  one  body  relative  to  the  other  (since 
it  is  relative  motion,  a  non-inertial,  non-rotating  coordinate  system  may  be  used  with  its 
origin  in  the  central  body)  [23:14]: 


d2r  a  _ 

+  ~r  =  ° 


(10) 


dt2  rz 

where  /v  is  a  gravitational  parameter  equivalent  to  G(mi+m2).  This  value  varies  based  on 
the  attracting  bodies.  When  the  mass  of  the  central  body  (mO  is  significantly  larger  than 
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its  satellite  (m2),  the  satellite’s  mass  can  be  ignored  and  /./  is  simply  G(mi).  For  earth  this 
value  is  approximated  as  3. 986012x10s  km3/s2  [23:429]. 

The  two-body  equation  is  a  useful  tool  for  a  quick  approximation  of  more  complex 
systems,  but  its  accuracy  is  limited  by  the  assumptions  used  for  its  derivation.  It  assumes 
the  bodies  are  perfectly  spherical  and  symmetric  in  composition  so  that  they  may  be 
treated  as  simple  point  masses  in  the  analysis  [23:11].  In  reality,  the  nonspherical  shape 
of  each  body  introduces  variations  in  gravitational  forces.  The  two-body  equation  further 
assumes  that  the  only  force  present  on  each  body  is  a  gravitational  force  acting  along  a 
line  joining  the  centers  of  each  body  [23:12],  This  discounts  the  effects  of  perturbations 
such  as  drag,  solar  radiation  pressure,  the  gravitational  influence  of  additional  bodies,  etc. 
Each  of  these  perturbations  must  be  accounted  for  if  a  higher  accuracy  solution  is  desired. 

2.1.2  Integrals  (Constants)  of  Motion 

The  solution  to  the  two  body  problem,  like  any  differential  equation,  requires 
constants  of  integration  known  as  integrals  of  motion,  or  more  simply,  just  integrals. 
Since  the  two-body  problem  consists  of  three  second  order  differential  equations  for  each 
body,  the  order  of  the  system  or  degrees  of  freedom  are  12  [given  by  (3x2)  x  2  bodies] 
[24:37].  To  reduce  the  order  to  zero,  12  integrals  are  needed. 

We  can  start  by  reducing  the  system  to  three  second-order  differential  equations  using 
conservation  of  linear  momentum  to  provide  the  first  six  constants  from  the  initial 
position  and  velocity  components  of  the  system’s  center  of  mass  [24:37].  This  transforms 
the  origin  to  the  barycenter  of  the  central  body  since  we  assume  the  satellite’s  mass  is 
negligible  [24:37].  The  remaining  six  are  used  to  define  the  shape  and  size  of  the 
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satellite’s  orbit  about  the  central  body,  the  orientation  of  the  orbit  using  Euler  angles,  and 
the  position  of  the  satellite  on  the  orbit  [25:208].  Since  angular  momentum  is  conserved, 

three  constants  come  from  the  components  of  the  angular  momentum  vector,  h  [25:80]. 
Three  more  constants  come  from  Kepler’s  first  law  which  gives  the  constant  eccentricity 

vector,  e.  Since  h  is  defined  perpendicular  to  e,  the  constraint  of  h-e  =  0  means  the 
components  of  angular  momentum  and  eccentricity  only  contribute  five  independent 
constants  of  integration  [25:208].  Kepler’s  second  law  can  be  used  with  the  polar 
component  of  the  angular  momentum  vector  to  provide  the  final  constant  of  integration  - 
the  time  of  perigee  passage  [24:37].  The  conservation  of  energy  also  provides  a  constant 
of  integration  [24:37],  but  it  is  not  independent  of  angular  momentum  and  eccentricity. 

Given  the  final  six  integrals,  the  satellite’s  orbit  can  be  completely  specified  in  space 
with  six  orbital  elements.  The  eccentricity  and  angular  momentum  scalar  quantities 
define  the  orbit  in  the  plane  [25:208].  The  inclination,  right  ascension  of  the  ascending 
node  and  argument  of  perigee  are  the  three  Euler  angles  that  define  the  orientation  of  the 
orbit.  The  first  two  are  directly  mapped  from  the  angular  momentum  vector  while  the 
latter  can  be  mapped  from  both  the  angular  momentum  vector  and  eccentricity  vector. 
Finally,  a  point  in  the  orbit  is  defined  by  the  true  anomaly  which  leads  to  the  time  since 
perigee  passage. 

For  //-body  problems,  n  >  3,  the  solution  requires  6 n  integrals  of  motion  [24:37].  Six 
are  given  by  the  conservation  of  linear  momentum,  one  is  given  by  the  conservation  of 
energy,  and  three  more  are  given  by  conservation  of  angular  momentum,  for  a  total  of  ten 
constants  [24:37].  The  additional  integrals  obtained  with  Kepler’s  laws  in  the  two-body 
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problem  are  not  applicable  for  n  >  3,  so  the  system  can  only  be  reduced  to  order  6n  -  1 0 
[24:37].  Thus  a  closed  form  solution  for  systems  greater  than  two  bodies  is  not  possible, 
but  simplifications  to  the  general  problem  have  yielded  analytical  solutions  such  as  the 
restricted  three-body  problem  (not  discussed  here). 


2.1.2.1  Specific  Angular  Momentum 

In  a  Keplerian  system,  the  only  force  acting  on  an  orbiting  satellite  is  the  gravitational 
force  directed  radially  toward  the  larger  central  body.  Without  other  tangential  forces  to 
alter  the  rotational  motion,  angular  momentum  is  conserved.  A  constant  called  the 
specific  angular  momentum  can  be  derived  independent  of  mass  by  first  cross 
multiplying  the  two-body  equation  with  the  position  vector  [24:24]: 


_  ^  jU  _ 

r  x  r  +  r  x  —r  =  0 


(11) 


The  second  term  vanishes  since  r  x  r  =  0  and  the  first  term  is  just  the  differential 
[24:24]: 


^-(rxr)  =  rxr  +  rxr  =  rxr 
dt v  J 


(12) 


Substituting  (12)  back  into  (11)  and  integrating  yields  a  vector  constant  of  integration, 
namely  constant  specific  angular  momentum  where  “specific”  signifies  the  mass 
independence  of  the  equation. 


h  =  r  x  v  = 


yvz  -  zvy 


[xvy  -  yvx\ 


=  constant 


(13) 


Since  h  is  a  constant  vector  perpendicular  to  r  and  v,  the  satellite’s  motion  is  fixed  in 
the  plane  containing  r  and  v  called  the  orbital  plane  [23:17], 
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2.1.2.2  Kepler’s  First  Law,  the  Trajectory  Equation,  and  Eccentricity 

Kepler’s  first  law  states  that  planets  follow  elliptical  paths  with  the  sun  at  one  focus 
[24:10].  This  statement  can  be  extended  to  include  all  conic  sections:  circles,  ellipses, 
parabolas  and  hyperbolas  [24:29].  The  mathematical  representation  of  Kepler’s  first  law 
describing  orbital  motion  on  a  conic  section  is  called  the  trajectory  equation. 

To  derive  the  trajectory  equation,  start  by  cross  multiplying  the  two-body  equation 
into  the  specific  angular  momentum  vector  [24:27]: 

r  x  ft  +  -^-T  x  ft  =  0  (14) 

y  J 

Since  angular  momentum  is  constant  ( ft  =  0),  the  first  term  can  be  written  [24:27]: 

^(r  x  ft)  =  r  x  ft  +  r  x  ft  =  r  x  ft  (15) 

Substituting  the  definition  of  angular  momentum,  (13),  into  the  second  term  in  (14) 
[24:28]: 

—  r  x  ft  =  —  |r  x  (r  x  r) J  (16) 

Using  the  “bac  -  cab”  rule  for  the  cross-products  [24:28]: 

^rxh  =  ^[r(f-r)-r(r-r)]  (17) 

Further  simplification  with  the  important  identities  r  ■  r  =  rr  and  r  ■  r  =  r2  yields 
[24:28]: 


i u  _  ^  /j.  r_ 

pj-r  x  ft  =  —  [ r(rr )  —  r(r2)J  =  — rr  —  —v 


(18) 


Recognizing  this  as  the  time  derivative  [24:28]: 
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(19) 


d  fr\  rr  —  rr  /j.  _  /i  _ 

— ju  — I  -  =  -jtf - ^ —  =  — rrr - 17 

dt  \r )  rt  rt  r 

Substituting  (15)  and  (19)  back  into  (14)  and  rearranging  [24:28]: 

4(?xft-„t)  =  0  (20) 

Integrating  [24:28]: 

_l  ^  r  ^ 

r  x  h  =  u-  +  B  (21) 

r 

B  is  the  vector  constant  of  integration,  sometimes  referred  to  as  the  Laplace  vector 
(named  after  French  mathematician  Pierre-Simon  Laplace)  [25:78],  which  is  directed 
toward  the  closest  point  to  the  central  body  called  periapsis. 

To  get  a  scalar  relation  for  (21),  dot  multiply  both  sides  by  r  [24:28]: 

r  ■  (r  x  h)  =r-^j.-  +  B^j  (22) 

The  left-side  is  simplified  using  the  vector  identity  A  ■  (B  x  C)  =  (A  XB)-C  from 
which  we  obtain: 


r  ■  (r  x  h)  =  (r  x  r)  ■  h  =  h-  h  =  h2  (23) 

Substituting  this  expression  and  expanding  the  right  side: 

h2  =  (r^l  +  r-B  (24) 

r 

From  the  definition  of  the  dot  product,  r  ■  r  =  r2  and  r  ■  B  =  rB  cos  (0)  where  0  is  the 
angle  between  the  fixed  vector  B  and  the  variable  position  vector  r.  Substituting: 

h2  =  /rr  +  rB  cos(0)  (25) 


Solving  for  position: 
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r  = 


(26) 


hi 

1  +  cos(6 ) 

This  expression  is  identical  to  the  polar  equation  of  a  conic  section: 

V 

r  = - 

1  +  e  cos(v ) 

Three  important  observations  between  equations  (26)  and  (27)  are  that  eccentricity,  e,  is 
equal  to  B/p,  the  parameter  or  semilatus  rectum,  p,  is  h2  /p,  and  the  angle  v,  called  true 
anomaly,  represents  the  angle  from  periapsis  to  the  current  location  on  the  orbit.  The 
expression  given  by  (27)  is  known  as  the  trajectory  equation  or  orbit  equation. 

Given  the  relationship  between  eccentricity  and  the  Laplace  vector  in  the  trajectory 
equation,  we  can  rewrite  (21)  as  the  eccentricity  vector  -  a  constant  of  integration: 


_  v  x  h  r 
e  = - 


(28) 


p  r 

Like  the  Laplace  vector,  the  eccentricity  vector  points  toward  periapsis  and  its  magnitude 
reveals  the  shape  of  the  Keplerian  orbit.  The  latter  property  will  be  discussed  more  in 
section  2.1.3. 


2.1.2.3  Kepler’s  Second  and  Third  Laws 

Kepler’s  second  law  states  that  the  line  joining  the  planet  to  the  sun  sweeps  out  equal 
areas  of  the  orbit  in  equal  time  [24:10],  The  mathematical  representation  of  this  over  a 
differential  element  of  area  along  the  conic  section  is  given  by  [24:30], 


dt  = 


dA 

h/2 


(29) 
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where  h! 2  is  the  areal  velocity.  Integrating  over  2n  radians  of  v  yields  the  period  of  one 


orbit  -  another  constant  of  motion: 


P  = 


2nab 

h 


(30) 


Notice  the  area  of  an  ellipse  is  given  by  nab  where  a  is  the  semimajor  axis  and  b  is  the 
semiminor  axis.  If  the  orbit  was  circular,  it  would  be  na2.  The  period  for  a  parabolic  or 
hyperbolic  orbit  must  be  derived  separately. 

It  is  much  more  insightful  to  alter  the  form  of  (30)  by  substituting  the  relation  for  the 
semiminor  axis,  b  =  aVl  —  e2  =  Jap,  where  the  semilatus  rectum,  p,  is  h2 /p: 


(31) 


This  is  precisely  the  expression  of  Kepler’s  third  law:  “The  square  of  the  period  of  a 
planet  is  proportional  to  the  cube  of  its  mean  distance  to  the  Sun  [24:10].” 

Kepler  used  the  previous  result  to  define  the  mean  angular  rate  of  a  planet  about  the 
Sun  in  units  of  radians  per  unit  time  as: 


2n  l~p~ 

n  P  a3 

The  constant  gravitational  parameter,  p,  can  be  solved  to  reveal  the  correlation: 

n2a3  =  p  =  Gim-t  +  m2) 


(32) 


(33) 


2.1.2.4  Specific  Mechanical  Energy  or  Vis-Viva  Equation 

A  satellite’s  total  mechanical  energy  -  the  sum  of  its  kinetic  and  potential  energy  -  is 
conserved  in  a  gravitational  field  void  of  external  energy  sources  and  sinks.  As  the 
satellite  orbits  a  central  body  it  continually  exchanges  one  form  of  energy  for  the  other.  A 
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constant  called  the  specific  mechanical  energy  can  be  derived  independent  of  mass  by 
first  dot  multiplying  the  two-body  equation  with  the  velocity  vector  [24:26]: 

v  ■  v  +  v  ■  -^-r  =  0  (34) 

y  J 

The  vectors  can  be  eliminated  by  recognizing  the  general  form  of  the  dot  product 
r  -v  =  rv  cos  9  and  realizing  that  the  radial  component  of  velocity,  not  to  be  confused 
with  the  magnitude  of  velocity,  is  r  =  v  cos  6  [24:26]: 


/r 

vv  +  r  —z  r  =  0 


Noticing  the  presence  of  the  following  time  derivatives  [24:26]: 

d  (v2\ 


dt  \  2 


=  vv 


d  /  ii\  _  j u 

dt  v  rf  r2 


(35) 


(36a) 

(36b) 


Substituting  into  (35)  yields  [24:26]: 


d 

dt 


(37) 


Integrating  both  sides  produces  the  energy  integral  or  vis-viva  (“living  force”)  equation 
[24:26]: 


v2  u 

£=— - +  C 

2  r 


(38) 


The  first  term  is  the  kinetic  energy.  The  second  term  is  potential  energy.  The  constant  c  is 
arbitrary  and  defined  by  the  physics  community  as  ii/R@,  but  in  astrodynamics  it  is 
defined  as  zero  [24:26]. 
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For  elliptical  orbits  that  include  rectilinear  ellipses,  another  form  of  the  energy 
equation  can  be  formed  in  terms  of  the  semimajor  axis,  a.  Start  by  dot  multiplying  the 
eccentricity  vector  by  itself  [26: 116]: 

e2  =  e  ■  e  =  —  (v  x  h)  ■  (v  x  /i)  —  — r  ■  (v  x  h)  +  1  (39) 


Simplification  yields  [26:116]: 


_  h2  (2  v2) 

1  —  e2  =  — ( - 


(40) 


M  \r  !i  J 

Recognizing  the  semilatus  rectum,  p  =  h2  /  jx,  and  knowing  from  the  geometry  of  an 
ellipse  that  p  =  a(l  —  e2),  the  second  term  on  the  right-hand  side  is  the  inverse  of  the 
semimajor  axis,  a  [26:116]: 


«  =  (-— -f 

\r  nj 

This  can  be  rearranged  to  match  the  right-hand  side  of  the  energy  equation  from  (38): 

ix  V2  IX 

2a  2  r 

Plugging  (42)  into  (38)  gives  the  alternate  form  of  the  energy  equation  [24:27], 

p 


(41) 


(42) 


£  =  — 


2  a 


(43) 


which  can  also  be  rearranged  to  describe  the  semimajor  axis: 

a  =  -TT 

2e 

The  traditional  form  of  the  vis-viva  equation  is  given  by  solving  (42)  for  v2  [24:27]: 

>2  In 


(44) 


v‘ 


=  /<(-—) 


(45) 
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2.1.3  Kepler’s  Equation 

The  final  constant  of  motion  mentioned  previously  as  the  time  of  periapsis  passage,  T, 
is  defined  once  every  orbital  period  when  the  satellite  passes  closest  to  the  central  body. 
Alone,  it  is  merely  a  point  in  time,  but  its  combination  with  the  current  time,  the  time 
since  periapsis  passage,  represented  by  t  —  T,  can  be  used  to  identify  the  satellite’s 
location  along  the  orbit.  Kepler  sought  a  method  that  relates  this  change  in  time  with  the 
angular  displacement  in  an  orbit.  His  second  law  states  that  equal  areas  are  swept  out  in 
equal  times,  so  he  resolved  that  all  time-to-area  ratios  can  be  equated  to  the  ratio  for  one 
complete  orbit  [24:51]: 


t-T 


(46) 


A1  nab 

The  unknown  area,  A;,  is  determined  from  geometry  by  encompassing  an  ellipse  with  an 
auxiliary  circle  shown  in  Figure  3.  First,  visually  define  the  area,  A;  [24:52]: 

At  =  Area  PCB  -  A2  (47) 

Triangle  A2  is  determined  by  substituting  trigonometric  relations  into  the  familiar 
geometric  area  of  a  triangle,  base  *  height /2  [24:52]: 

1  (b  \  ab 

A2  =  ~  (ae  —  a  cos(ii))  ( —  a  sin(F)  )  =  —  (e  sin (£)  —  cos(fi)  sin (£))  (48) 

2  \a  /  2 


The  area  of  the  ellipse  inside  points  PCB  is  determined  by  scaling  down  the  similar  area 
of  the  auxiliary  circle  inside  points  PCB’  by  b/a  (the  relationship  between  vertical 
distance  on  an  ellipse  and  a  circle  -  see  [25:183]).  The  area  of  PCB’  is  merely  the  circular 
area  of  POB’  minus  the  triangular  area  OB’C. 
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Figure  3.  Ellipse  circumscribed  in  a  circle  to  provide  geometry  for  Kepler’s  Equation. 


Together,  this  yields  [24:53]: 

b  ab 

-  =  —[£■  —  cos  (E)  sin(E)]  (49) 
a  2 


i4rea  PCB  = 


azE  1 

— - -(a  cos(£,))(asin(E)) 


Substituting  (48)  and  (49)  into  (47)  [24:53]: 

ab 

A1=—[E  -  e  sin(E)] 

Plugging  this  result  into  (46)  and  rearranging  gives  [24:53]: 

p  2n(t  -  T) 

E  —  e  sin(E) 


(50) 


(51) 
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Equating  the  two  expressions  for  the  period  of  a  satellite  given  by  (31)  and  (51)  yields 
Kepler ’s  equation  -  a  powerful  tool  for  determining  time  since  periapsis  passage  from  the 
eccentric  anomaly,  semimajor  axis  and  eccentricity  [24:53]: 


M 


aJ 


(t-T) 


(52) 


p  E  —  e  sin  (£) 

Kepler  rearranged  (52)  to  define  a  parameter  called  the  mean  anomaly,  M,  which 
represents  the  theoretical  angle  swept  out  by  a  satellite  travelling  at  its  mean  angular 
velocity  on  a  circular  orbit  of  radius  a  [24:53]: 


M  =  E  -  e sin(E)  =  4( t-T ) 

vJ  a6 


(53) 


The  mean  anomaly  has  no  physical  significance,  but  its  importance  comes  from 
relating  the  eccentric  anomaly  with  time  [24:54].  First  calculate  mean  anomaly  from  time 
and  the  semimajor  axis.  Then  solve  the  eccentric  anomaly  from  the  mean  anomaly’s 
transcendental  form.  Once  E  is  solved  using  graphical  or  numerical  methods,  it  can  be 
related  to  its  corresponding  true  anomaly.  Start  by  defining  the  eccentric  anomaly  as  a 
function  of  true  anomaly  using  the  auxiliary  circle  from  Figure  3  [24:55]: 

rsin(u) 


sin  (£)  = 


cos(E)  = 


ae  +  r  cos(u) 


a 


(54) 

(55) 


The  trajectory  equation  (27)  can  be  substituted  for  the  position  r  to  give  [24:55]: 

a(l  —  e2)  sin(u)  sin(u)  Vl  —  e2 


sin(E)  = 


a(l  +  e  cos(u))Vl  —  e2  1  +  e  cos(u) 


(56) 
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cos  (£)  =  e  + 


g( 1  —  e2) 

(1  +  e  cos(u)) 
a 


cos(u) 


e  +  cos(u) 
1  +  e  cos(u) 


(57) 


Solving  (57)  for  cos(u)  [24:55]: 


cos(u)  = 


cos(£)  —  e 
1  —  e  cos  (E) 


(58) 


Next  determine  an  expression  for  position  as  a  function  of  semimajor  axis  and 


eccentricity  by  rearranging  (55)  [24:55]: 

acos(E)  —  ae 

7"  =  - 

cos(u) 


(59) 


Plug  (58)  into  (59)  to  eliminate  true  anomaly  [24:56]: 


a  cos  (£)  —  ae 

cos  (£)  —  e 
1  —  e  cos  (£) 


=  a(l  —  e  cos(ij)) 


(60) 


Solve  (54)  for  sin(u)  [24:56]: 


a  sin(£’)  Vl  —  e2 
sin(u)  = - 


(61) 


Plug  (60)  into  (61)  to  solve  true  anomaly  in  terms  of  e  and  E  [24:56]: 


sin(u)  = 


sin(ii)  Vl  —  e2 
1  —  e  cos  (E) 


(62) 


Now  (58)  and  (62)  can  be  substituted  into  the  half-angle  identity  for  tangents  and 
simplified  as  [24:56]: 


sin(E)  Vl  —  e2 

_  1  -  e  cos(E)  _  sin(ff)  Vl  -  e2 

an  W  cos (£~)  —  e  cos (E)  —  e  +  1  —  e  cos (E) 

1  —  e  cos (£) 


(63) 


The  denominator  in  (63)  is  simply  the  product  of  (1  —  e)(cos (E)  +  1)  [24:56]: 
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(64) 


sin(ii)  Vl  —  e2 

tan  V2 )  =  (1  —  e)(cos(£’)  +  1) 

Again  the  tangent  half  angle  can  be  used  for  the  eccentric  anomaly  components  to 
simplify  (64)  [24:56]: 


tan 


©- 


1  +  e  ( E 

tan 


1  —  e 


(£) 


(65) 


Upon  rearrangement,  the  eccentric  anomaly  is  then  just: 


E  =  2  tan  1 


1  —  e 

- - tan .  „ 

1  +  e  V2 


© 


(66) 


2.1.4  State  Representation 

The  state  of  a  satellite  can  be  completely  specified  in  time  by  just  six  quantities.  In  a 
standard  Cartesian  reference  frame,  these  are  given  by  the  position  and  velocity  vectors. 
Additionally,  an  element  set  can  be  used  that  consists  of  scalar  magnitudes  and  angles. 
The  most  commonly  used  one  is  referred  to  as  classical  or  Keplarian  orbital  elements.  Its 
canonical  counterpart  is  known  as  the  Delaunay  elements.  Although  the  classical  orbital 
elements  (COEs)  are  more  intuitive  to  envision  than  a  Cartesian  set,  they  do  suffer  from 
some  orbital  geometry  challenges.  To  bypass  these  difficulties,  variations  of  the  COEs 
exist,  such  as  the  equinoctial  elements  and  their  canonical  counterparts  known  as 
Poincare  elements.  Nevertheless,  the  complete  bill  of  fare  is  beyond  the  purview  of  this 
study,  so  we  restrict  ourselves  to  a  review  of  the  COEs  and  their  Delaunay  elements. 
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2.1.4.1  Classical  Orbital  Elements 


The  state  vector  typically  used  to  identify  a  two-body  orbit  in  space  is  given  by  a  set 
of  six  elements  called  the  COEs  or  Keplarian  elements  which  include  the  semi  major 
axis,  a,  the  eccentricity,  e,  the  inclination,  i,  argument  of  perigee,  co,  right  ascension  of 
the  ascending  node,  Q,  and  true  anomaly,  v  [24:104],  The  COEs  are  the  initial  conditions 
required  to  solve  the  initial  value  problem  for  two-body  orbital  motion  [24:103],  The 
COEs  can  be  formulated  from  the  position  and  velocity  vectors  in  rectangular  coordinates 
typically  obtained  from  observation  or  measurement. 


perigee  direction 


sale Hue  ’s  position 


X 

/ 

Y 

ve  rnal  equinox 


Figure  4.  Classical  orbital  elements.  Credit:  NASA 


The  first  element  -  the  semimajor  axis,  a  -  defines  the  size  of  the  orbit  [24:104], 
Some  texts  supplant  this  element  with  specific  angular  momentum,  h,  or  the  semilatus 
rectum,  p,  since  the  semimajor  axis  is  infinite  for  parabolic  orbits  [25:209],  Derived  from 
the  semimajor  axis,  the  mean  motion,  n,  can  also  be  used  to  define  the  shape  of  any  orbit 
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except  for  parabolic  orbits  [24:105].  For  non-parabolic  orbits,  the  semimajor  axis  is  most 


easily  calculated  from  the  energy  equation  or  more  expressly  (41). 

The  second  element  -  the  eccentricity,  e  -  defines  the  shape  of  the  orbit  and  is  the 
magnitude  of  the  eccentricity  vector  [24:105].  Start  by  expanding  the  vector  triple  cross 
product  from  (28)  to  get  [24:106]: 

_  (v  •  v)r  —  (r  •  v)v  r 

e  = -  (67) 

/r  r 

Simplifying  the  dot  product  of  the  velocity  vectors  and  combining  position  vector  terms 
yields  the  final  equation  [24:106]: 

_  =  v)v  (68) 

The  magnitude  of  the  eccentricity  vector  is  simply: 


|e||  =  \ef  +  ef  +  el 


(69) 


For  non-parabolic  orbits  the  eccentricity  can  also  be  determined  by  [24:106]: 


e  = 


(70) 


M  a 

When  the  magnitude  of  eccentricity  is  zero,  the  orbit  is  circular.  For  0  <  e  <  1,  the 
orbit  is  elliptic.  When  its  magnitude  is  one,  the  orbit  is  parabolic  and  for  e  >  1,  the  orbit 
is  hyperbolic. 

The  third  element  -  inclination,  i  -  defines  the  tilt  of  the  orbit  plane  or,  more 

specifically,  the  angle  between  the  inertial  K  axis  and  the  angular  momentum  vector,  h 
[24:107],  It  is  also  helpful  to  think  of  it  as  the  angle  between  the  orbit  plane  and  the 
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equatorial  plane  ( 7/  plane)  of  the  central  body.  The  inclination  is  defined  without 
quadrant  ambiguity  (since  i  is  always  in  the  range  of  0°  to  180°)  as  [24:107], 


i  =  cos 


(71) 


where  h  is  given  by  (13). 

The  fourth  element  -  right  ascension  of  the  ascending  node  (RAAN),  Q  -  defines  the 
angle  from  the  inertial  7  axis  to  the  point  on  the  equatorial  plane  where  the  satellite 
crosses  from  the  southern  hemisphere  into  the  northern  hemisphere  [24:107].  This  point 
is  called  the  ascending  node  and  it  is  located  in  the  direction  of  the  nodal  vector  given  by 
[24:107]: 


^  _  \~hj] 

n  =  Kxh  =  ^ 

-  0  - 


(72) 


The  RAAN  is  then  defined  as  [24:107]: 


fl  =  cos  1 


(73) 


A  quadrant  check  is  necessitated  for  (73).  The  RAAN  will  range  from  0°  <  fl  <  180°  if 
the  ascending  node  is  on  the  positive  side  of  the  IK  plane,  meaning  the  j  component  of 
the  nodal  vector  is  greater  than  zero,  nj  >  0.  If  the  ascending  node  is  on  the  negative  side 
of  the  IK  plane,  meaning  iij  <  0,  then  the  correct  angle  and  quadrant  are  given  by 
D.  =  360°  —  fl  so  that  RAAN  will  range  from  180°  <  fl  <  360°  [24:107]. 

We  could  have  also  defined  RAAN  as: 


32 


n  =  sin  1 


Jn  \ 

\j\\\\n\\) 


(74) 


Combining  (73)  and  (74),  we  can  define  RAAN  with  an  ATAN2  function  to  avoid 
quadrant  ambiguity  altogether: 


n  =  tan"1(|)  =  tan"1(-|)  <75) 

The  fifth  element  -  argument  of  perigee,  co  -  measures  the  angle  from  the  ascending 

node  to  periapsis.  Since  the  eccentricity  vector  points  to  periapsis,  the  argument  of 
perigee  is  simply: 


co  =  cos  1 


n  ■  e 

PiiPii, 


/ 


=  cos  1 


htej  -  hjej 

\e  lhf  +  hf 
\  a/  J ' 


(76) 


A  quadrant  check  is  needed  since  co  can  range  from  0°  to  360°.  If  periapsis  lies  below  the 
equatorial  plane,  i.e.  the  “k”  component  of  the  eccentricity  vector  is  less  than  zero,  ty  <  0, 
the  correct  angle  and  quadrant  are  given  by  a)  =  360°  —  co  [24:108]. 

The  sixth  element  -  true  anomaly,  u  -  measures  the  angle  between  periapsis  and  the 
position  vector  which  indicates  the  current  location  of  the  satellite  in  the  orbit,  so: 


o  =  cos  1 


e  r 

MM, 


(77) 


A  quadrant  check  is  needed  since  v  can  range  from  0°  to  360°.  If  the  satellite  is  moving 
toward  periapsis  (away  from  apoapsis),  indicated  when  r  ■  v  <  0,  then  the  correct  angle 
and  quadrant  are  given  by  u  =  360°  —  u  [24:108], 
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Degeneracies  result  from  the  standard  six  elements  when  the  orbit  is  equatorial, 
circular  or  both.  Under  these  geometric  conditions,  alternate  elements  must  be  used  to 
define  the  orbit. 

For  elliptical  equatorial  orbits  where  i  =  0°  or  180°,  the  nodal  vector  is  undefined 
causing  the  RAAN  to  exhibit  a  singularity.  The  solution  is  to  replace  the  RAAN  with  the 
true  longitude  of  periapsis,  G)true  ,  which  is  the  angle  measured  positive  eastward  from 
the  vernal  equinox  to  the  eccentricity  vector  [24:109].  Mathematically  this  is  defined  as: 


to true  COS 


(78) 


A  quadrant  check  is  necessitated  since  (otrue  values  range  from  180°  to  360°  (or  similarly 
-180°  to  0°).  If  the  “j”  component  of  the  eccentricity  vector  is  less  than  zero,  ej  <  0,  then 
the  correct  angle  and  quadrant  are  given  by  6otrue  =  360°  —  ojtrue. 

A  note  is  warranted  on  the  difference  between  true  longitude  of  periapsis,  cotrue,  and 
a  similar  quantity  called  longitude  of  periapsis,  to,  given  by  [24:109]: 


to  —  fl  +  to 


(79) 


Longitude  of  periapsis  defines  the  angle  from  the  primary  axis  to  periapsis  by  summing 
two  angles  in  different  planes,  whereas  true  longitude  of  periapsis  is  a  measure  of  the 
angle  between  the  primary  axis  and  the  eccentricity  vector  [24:109].  At  small  inclinations 
they  are  similar,  but  as  inclination  increases  the  difference  grows  [24:110], 

For  circular  inclined  orbits,  the  eccentricity  vector  is  undefined  and  a  singularity 
results  when  trying  to  determine  argument  of  perigee  and  true  anomaly.  The  angle  used 
instead  is  the  argument  of  latitude,  u,  which  is  measured  positive  from  the  ascending 
node  to  the  satellite’s  position: 
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It  =  cos 


(80) 


/  n  r  \  _  ^(Ky-hyxX 

(||H||||F||)“C0S  \rJi% Th}J 

A  quadrant  check  is  needed  since  u  can  range  from  0°  to  360°.  If  the  satellite  is  south  of 
the  equatorial  plane,  indicated  when  rk  <  0,  then  the  correct  angle  and  quadrant  are  given 
by  it  =  360°  -u  [24:108], 

When  argument  of  perigee  and  true  anomaly  are  defined  for  non-degenerate  cases, 
the  argument  of  latitude  is  just: 


u  =  o)  +  v  (81) 

The  last  degenerate  case  occurs  when  an  orbit  is  circular  and  equatorial  [24:111],  In 
this  case,  the  ascending  node  vector  and  eccentricity  vector  are  undefined  causing 
singularities  in  RAAN,  argument  of  perigee  and  true  anomaly.  The  angle  used  instead  is 
the  true  longitude,  Atrue,  measured  eastward  from  the  primary  axis  to  the  satellite’s 
current  position  [24:  111]: 


A true  COS 


(82) 


A  quadrant  check  is  needed  since  Atrue  can  range  from  0°  to  360°.  If  the  “j”  component 
of  the  position  vector  is  less  than  zero,  r;  <  0,  then  the  correct  angle  and  quadrant  are 
given  by  Atrue  =  360°  -  Atrue  [24:111], 

For  non-degenerate  cases,  the  true  longitude  can  be  approximated  as  [24:  111]: 

A true  ~  H  4"  (O  +  V  (83) 

It  is  only  approximate  since  RAAN  is  in  a  different  plane  than  argument  of  perigee  and 
true  anomaly  in  an  inclined  orbit  [24:  111], 
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2.1.4.2  Delaunay  Elements 

The  classical  orbital  elements  have  the  benefit  of  being  a  relatively  intuitive  way  of 
viewing  an  orbit,  but  they  make  the  equations  of  motion  more  difficult  to  solve  than 
another  choice  of  variables.  Their  canonical  equivalent — one  which  satisfies  Hamilton’s 
equations  (discussed  in  §2.2.2) — renders  more  trivial  equations  of  motion  with  its  use  of 
generalized  coordinates  and  conjugate  momenta.  For  their  simplicity,  these  so-called 
Delaunay  elements  are  routinely  used  in  perturbation  theory.  Their  derivation  is  left  to 
reference  textbooks  such  as  Wiesel  [27:62];  however,  the  result  is  given  here  in  which  the 
capitals  denote  the  momenta: 


L  =  yfjia 

(84) 

G  =  LyJ  1  —  e2 

(85) 

H  =  G  cos  i 

(86) 

l  =  M 

(87) 

g  =  co 

(88) 

h  =  a 

(89) 

It  should  be  noted  that  the  Delaunay  elements  suffer  from  the  same  degeneracies 
when  the  eccentricity  and  inclination  approach  zero. 

2.1.4.3  Coordinate  Systems 

The  differential  equations  of  motion  in  both  Newtonian  and  Hamiltonian  systems 
require  an  inertial  celestial  reference  frame  for  their  solution.  The  closest  realization  of 
this  in  our  neighborhood  of  the  universe  is  the  extra-galactic  coordinate  system  anchored 
at  the  center  of  the  Milky  Way  Galaxy.  Unfortunately,  it  is  not  perfectly  inertial,  nor  is  it 
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realistic  for  Earth  orbiting  satellites.  A  sufficiently  inertial  reference  frame  can  be  used 
that  places  the  origin  at  the  Earth’s  center  at  any  epoch  of  interest,  directs  the  z-axis  along 
the  Earth’s  axis  of  rotation  in  the  right-handed  sense,  and  directs  the  x-axis  toward  the 
vernal  equinox  in  the  equator  of  epoch.  This  is  known  as  a  true  of  date  (TOD)  frame  as 
shown  in  Figure  5.  It  changes  very  slowly  over  time  due  to  planetary  precession,  luni- 
solar  precession,  and  nutation  from  the  Sun  and  the  moon;  however,  the  collective 
precessional  period  relative  to  the  fixed  stars  takes  nearly  26,000  years  [24:207].  As  such, 
the  apparent  diurnal  motion  of  the  vernal  equinox  is  very  close  to  the  apparent  motion  of 
the  fixed  stars.  This  makes  the  vernal  equinox  a  good  reference  point  for  a  sufficiently 
inertial  frame. 


EARTH’S  TRUE-OF-DATE 
ROTATIONAL  AXIS 


TRUE  EC 
D/ 


Figure  5.  True  of  Date  (TOD)  Cartesian  coordinate  system.  Credit:  NASA 
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Earth  orientation  models  are  used  to  determine  the  true  of  date  equator  (current 
position  of  the  equatorial  plane)  and  the  true  equinox  of  date  (intersection  of  the  current 
equatorial  plane  and  ecliptic).  The  latest  of  these,  published  in  2003,  is  known  as  IAU 
2000,  but  the  classical  Fundamental  Katalog  5  (FK5)  theory  will  be  exclusively  covered 
here  since  it  is  used  for  the  ISS — the  source  of  data  for  this  study. 

2.1.4.3.1  Pseudo-Inertial  J2000  Frame 

The  TOD  frame  just  introduced  is  not  a  good  tool  for  comparing  data  over  a  long 
period  of  time  since  the  frame  is  constantly  changing  (due  to  the  Earth’s  precession  and 
nutation),  albeit  very  slightly.  The  solution  for  high  precision  comparisons  is  to  transform 
all  data  into  a  TOD  frame  at  a  reference  epoch.  The  J2000  coordinate  frame  in  FK5 
theory  does  this  by  establishing  a  standardized  Earth-centered  inertial  (ECI)  frame  at 
epoch  January  1,  2000,  noon  TT.  It  places  the  origin  at  the  Earth’s  center  with  the  xy- 
plane  in  the  Earth’s  true  of  date  (TOD)  equator  at  epoch  J2000,  the  x-axis  directed 
toward  the  true  of  date  equinox  at  epoch  J2000,  and  the  z-axis  pointed  North  along  the 
Earth’s  rotational  axis. 

2.1.4.3.2  True  of  Date  Rotating  (Earth-Fixed  Greenwich)  Frame 

In  addition  to  celestial-based  inertial  reference  frames,  observations  of  satellites 
require  terrestrial-based  coordinate  frames  fixed  to  the  rotating  Earth.  These  come  in  the 
form  of  geocentric,  barycentric  and  topocentric,  but  only  geocentric  is  covered  here.  The 
geocentric  version  is  routinely  called  the  Earth-centered,  Earth-fixed  (ECEF)  frame 
which  means  that  its  orientation  remains  unchanged  with  respect  to  the  Earth’s  crust.  A 
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strict  definition  of  the  ECEF  frame  must  account  for  polar  motion — the  difference 
between  the  present  axis  of  rotation  called  the  Celestial  Ephemeris  Pole  (CEP)  and  the 
International  Reference  Pole  (IRP)  which  was  agreed  upon  in  1900  and  1905  [24:207]. 
The  resulting  coordinate  system  is  called  the  International  Terrestrial  Reference  Frame 
(ITRF),  which  U.S.  Air  Force  Space  Command  calls  the  Earth-centered  rotating  (ECR) 
frame  [24:158].  While  the  ITRF/ECR  is  great  for  surveys  or  navigation,  it  is  not 
conducive  for  satellite  observations  since  it  is  not  referenced  to  the  CEP. 

The  ideal  frame  for  observations  which  is  referenced  to  the  CEP  is  simply  a  rotating 
version  of  the  TOD  frame.  NASA  calls  this  the  true  of  date  rotating  (TDR)  frame, 
whereas  the  U.S.  Air  Force  Space  Command  calls  it  the  Earth-fixed  Greenwich  (EFG) 
frame  [24:158].  Just  like  the  TOD  frame,  the  origin  of  the  TDR  frame  is  fixed  to  the 
center  of  the  Earth,  but  the  x-axis  of  the  TDR  frame  is  fixed  to  a  local  meridian  at  the  true 
equator  of  date  instead  of  pointing  toward  the  true  equinox  of  date.  Generally,  the  prime 
meridian  is  used  and  the  z-axis  is  the  same  as  the  TOD  frame.  See  Figure  6  for  its 
depiction.  The  only  difference  between  the  TDR  frame  and  the  ITRF/ECR  frame  is  the 
correction  for  polar  motion  which  is  very  small  since  the  maximum  variations  between 
the  CEP  and  IRP  are  approximately  9  meters  [24:208]. 

Converting  from  TOD  to  TDR  and  viceversa  requires  knowledge  of  Earth’s  rotation 
relative  to  the  vernal  equinox — a  concept  known  as  sidereal  time.  Viewed  from  the  North 
Pole,  the  Local  Sidereal  Time,  0LST,  is  the  angle  measured  counterclockwise  from  the 
vernal  equinox  (defined  at  the  equator)  to  the  local  meridian.  When  the  prime  meridian 
(0°  longitude)  at  Greenwich  is  used,  there  are  two  forms  of  the  sidereal  time:  Greenwich 
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Figure  6.  True  of  Date  Rotating  (TDR)  Cartesian  coordinate  system.  Credit:  NASA 


mean  sidereal  time,  dGMST,  and  Greenwich  apparent  sidereal  time,  0GAST.  The  mean 
sidereal  time  is  measured  along  the  true  equator  from  the  mean  equinox  to  the  Greenwich 
meridian,  whereas  the  apparent  sidereal  time  is  measured  along  the  true  equator  from  the 
true  equinox  to  the  Greenwich  meridian.  In  other  words,  the  apparent  sidereal  time 
includes  the  secular  precession  of  the  equinox  and  the  periodic  nutation  effects,  whereas 
the  mean  sidereal  time  only  includes  the  secular  precession  of  the  equinox.  This 
difference  amounts  to  the  equation  of  the  equinoxes  [24:223], 

EQ equinox  =  AH'  cos(<T)  +  0.00264"  sin(ftmoon)  +  0.000063  sin(2ftmoon)  (90) 
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where  AT*  is  the  nutation  in  longitude,  e  is  the  mean  obliquity  of  the  ecliptic,  and  f lmoon 
is  the  RAAN  of  the  mean  lunar  orbit.  The  general  expression  for  Greenwich  apparent 
sidereal  time  is  then: 


@GAST  @GMST  equinox  (91) 

There  are  a  number  of  ways  for  calculating  the  Greewich  mean  sidereal  time,  but  a 
convenient  method  for  computers  is  given  by  [24:191], 


eGMST  =  67310.54841s  +  (876600^  +  8640184.812866S)7,[/T1 


(92) 


+  0.0931047^  -  6.2  X  1(T67^t1 
where  Tuti  is  the  number  of  Julian  centuries  since  epoch  J2000. 

Since  the  contribution  from  the  equation  of  the  equinoxes  is  on  the  order  of  a 
thousandth  of  a  degree,  the  apparent  sidereal  time  is  approximately  equal  to  the  mean 
sidereal  time.  For  many  applications,  the  mean  sidereal  time  is  close  enough;  however, 
this  author  will  incur  the  extra  complications  to  have  that  precision. 

Given  the  apparent  sidereal  time,  a  simple  rotation  about  the  z-axis  is  required  to 
convert  from  TOD  to  TDR: 


rTDR  ~  ^3  [^GASt\TTOD 
VTdR  =  ^3WgASt]VTOD  ~  X  rTDR 


(93) 

(94) 


2.1.5  Orbit  Perturbations 

In  the  Earth  system,  the  primary  perturbing  forces  include  gradients  in  the  Earth’s 
gravity  field,  atmospheric  drag,  solar  radiation,  and  gravitational  effects  from  the  sun  and 
the  moon.  Atmospheric  drag  (a  nonconservative  force)  decreases  exponentially  with 
altitude,  yet  dominates  at  altitudes  below  -315  km,  while  third-body  effects 


41 


(conservative)  and  solar  radiation  (nonconservative)  dominate  at  altitudes  above  -19,950 
km  [22:271].  In  between  these  altitudes,  the  gravity  gradient  (conservative)  is  the  largest 
of  the  perturbations.  It  is  this  region  upon  which  the  author  will  primarily  focus  since  the 
ISS  is  maintained  at  altitudes  near  350  km.  As  such,  third-body  effects  and  solar  radiation 
will  not  be  discussed  here. 


2.1.5.1  Nonspherical  Earth 

Kepler’s  equations  were  derived  with  the  assumption  that  the  central  body  is  a 
Newtonian  point  mass.  For  the  Earth  to  have  this  property  it  would  have  to  be  a  perfect 
sphere  with  uniform  mass  distribution  or  the  radius  would  have  to  be  large  enough  for  the 
Keplarian  term  to  dominate.  In  reality,  the  Earth  is  bulged  at  the  equator,  has  a  slight  pear 
shape,  and  is  flattened  at  the  poles  [28:142],  Geographical  features  also  contribute  to 
anomalies  in  the  gravitational  field  compared  to  a  uniform  featureless  surface.  Large 
concentrations  of  mass  such  as  mountain  ranges  add  to  the  gravity  field  whereas  ocean 
trenches  and  depressed  landmasses  reduce  the  gravity  field. 

Even  with  all  of  these  complexities,  the  Earth’s  geopotential  can  be  modeled  with 
extraordinary  accuracy  and  precision  with  a  spherical  harmonic  expansion  given  by 
[27:108]: 


w  /  \~n 

p"(cose) 


71  —  0  772  =  0 


(95) 


x  \Cnm  cos (m0)  +  Snm  sin(m0)] 

ju  is  the  Earth’s  gravitational  constant,  r  is  the  radius  from  the  Earth’s  center  to  the 
satellite,  R@  is  the  equatorial  radius  of  the  Earth,  n  is  the  degree  of  the  expansion,  m  is  the 
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order  of  the  expansion,  P™  are  the  Legendre  polynomials,  Cmn  and  Snm  are 
dimensionless  geopotential  coefficients,  6  is  the  colatitude  and  <p  is  the  geocentric 
longitude. 


Gravity  Anomaly  (mGal) 

Figure  7.  GRACE  gravity  map  depicting  gravity  anomalies  on  Earth.  Credit:  University 
of  Texas  Center  for  Space  Research  and  NASA 

The  coefficients  can  be  divided  into  three  classes.  The  coefficients  for  degree  m  =  0  are 
called  zonal  harmonics  and  depend  only  on  latitude  [27:113],  The  coefficients  for  which 
n  =  m  are  called  sectoral  harmonics  and  depend  only  on  longitude  [27:114],  Finally,  the 
coefficients  remaining  when  m  =£  0  and  n  =£  m  are  called  tesseral  harmonics  and  depend 
on  both  longitude  and  latitude  [27:1 15].  The  last  form  appears  as  a  “checkerboard  pattern 
of  regions  that  alternatively  add  to  and  subtract  from  the  two-body  potential”  [28:143]. 
The  coefficients  are  determined  experimentally  and  placed  in  a  gravity  model  matrix  such 
as  NASA’s  Earth  Gravity  Model  1996  (EGM96)  or,  more  recently,  models  produced 
from  the  Gravity  Recovery  and  Climate  Experiment  (GRACE). 


43 


2.1.5.2  Atmospheric  Drag 


The  basic  equation  for  atmospheric  drag  is  given  as  a  specific  force  or  acceleration 
[24:525], 


1  CjjA  2  vrei 
“drag  ~  ~2~^pVrel\^\ 


(96) 


where  cD  is  the  drag  coefficient,  A  is  the  cross-sectional  area  perpendicular  to  the  velocity 
vector,  m  is  the  satellite’s  mass,  p  is  the  atmospheric  density,  and  vrel  is  the  velocity 
vector  relative  to  the  rotating  atmosphere,  or  simply  [24:526]: 


_  dr  _  _ 

Vrel  =  di~M®Xr  = 


dx 


dy 


dz 


(97) 


When  working  in  the  hypersonic  realm  it  is  generally  common  to  characterize  the 
susceptibility  to  drag  with  the  so-called  ballistic  coefficient  (BC)  which  is  simply  an 
inverse  of  the  second  fraction  in  (96)  [24:525]: 


BC  = 


m 

cdA 


(98) 


A  satellite  with  a  larger  surface  area  and  larger  drag  coefficient  will  have  the  effect  of 
decreasing  the  BC  and  increasing  the  overall  drag.  Thus,  the  BC  has  an  inverse 
relationship  with  drag. 

Drag  is,  perhaps,  the  most  difficult  perturbation  to  model  because  it  necessitates 
knowledge  of  precise  attitudes  for  area  calculations  and  requires  a  good  understanding  of 
the  atmospheric  density  which  can  fluctuate  wildly  due  to  solar  flux,  geomagnetic 
variations,  and  the  molecular  structure  of  the  atmosphere  [24:526].  Changes  in  extreme 
ultraviolet  radiation  (EUV)  from  incident  solar  flux  cause  instantaneous  heating  of  the 
atmosphere.  This,  in  turn,  causes  density  changes  at  altitudes  above  the  level  heated  due 
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to  increased  particle  collisions.  Geomagnetic  storms  have  the  same  impact,  but  exhibit 
delayed  heating  with  density  variations  persisting  for  8  to  24  hours  after  ground-based 
magnetometers  indicate  the  storm  has  ended  [29:137]. 

Tascione  shows  that  for  a  minor  geomagnetic  storm  identified  by  the  planetary  Kp 
index  of  3,  a  satellite’s  in-track  displacement  (compared  to  predictions)  grew  to  5  nautical 
miles  in  just  over  12  hours.  For  a  severe  storm  identified  by  Kp  =  8,  a  satellite’s  in-track 
displacement  grew  to  50  nautical  miles  over  the  same  time  period.  Of  course,  polar 
orbiting  satellites  will  experience  the  most  severe  orbit  degradation  due  to  more  prevalent 
atmospheric  heating  near  the  poles,  but  density  changes  can  pervade  deep  into  lower 
latitudes  and  influence  all  satellite  orbits.  Skylab  was  a  perfect  example  of  this.  Even  with 
a  50  degree  inclination,  its  orbit  quickly  deteriorated  due  to  increased  solar  activity  at  the 
beginning  of  Solar  Cycle  21. 


2.1.6  Variation  of  Parameters  (VOP) 

The  six  classical  orbital  elements  (a,  e,  i,  ft,  a>,  T0)  were  shown  to  be  constants  of 
motion  in  a  Keplarian  system.  With  the  addition  of  time  as  a  seventh  parameter,  a 
satellite’s  position  and  velocity  could  be  precisely  determined.  In  reality,  all  of  these 
constants  are  changing,  some  faster  than  others.  This  is  due  to  perturbations  from  forces 
other  than  the  central  body’s  gravity.  These  additional  forces  alter  the  two-body 
equations  of  motion  given  by  (10)  to  yield  [24:578]: 


d2r  n  _  _ 

It2  =~^3r  +  apert 


(99) 
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One  method  for  solving  the  perturbed  orbit  is  called  variation  of  parameters.  VOP 
makes  the  assumption  that  the  perturbations  are  small  enough  that  the  solution  to  the 
unperturbed  two-body  system  will  suffice  for  the  perturbed  system.  Another  assumption 
is  that  the  constants  of  motion  from  the  unperturbed  system  can  be  represented  as  time 
varying  in  the  perturbed  system.  So  instead  of  integrating  the  rectangular  coordinates  in 
(99)  as  would  be  done  for  the  special  perturbation  technique  (Cowell’s  method),  VOP 
consists  of  six  first-order  differential  equations  describing  the  time  rate  of  change  of  the 
osculating  elements  (a,  e,  i,LL,  a>,  M).  The  term  osculating  is  used  because  the  elements 
are  no  longer  constant.  They  represent  the  satellite’s  position  at  one  instant  in  time  since 
they  are  continuously  altered  by  perturbations.  Nevertheless,  once  the  osculating 
elements  are  solved  via  the  differential  equations,  they  may  be  resolved  into  a  position 
and  velocity  vector  using  the  two-body  problem  solution.  This  process  repeats  itself  for 
each  point  in  time. 

There  are  two  well-known  techniques  for  VOP;  one  was  developed  by  Lagrange  and 
the  other  by  Carl  Friedrich  Gauss.  Lagrange’s  version  applies  to  conservative  forces 
while  Gauss’  also  works  for  non-conservative  forces  [24:577].  Both  are  used  as  a 
foundation  for  general  perturbation  theory  which  is  of  limited  interest  to  the  present 
study.  However,  Lagrange’s  version  will  be  developed  here  to  derive  Lagrange’s 
Planetary  Equations  of  Motion  which  contribute  a  set  of  secular  rates  to  the  basis 
frequencies  needed  for  spectral  analysis. 
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2.1.6.1  Lagrange  Planetary  Equations  -  Potential  Form 

Two  different  forms  of  Lagrange’s  Planetary  equations  can  be  derived:  one  arises  in 
an  acceleration  component  form  and  the  other  in  a  potential  form.  Since  the  author’s 
motivation  is  merely  to  come  up  with  a  set  of  secular  rates  from  the  J2  disturbing 
function,  only  the  potential  form  will  be  shown  here.  The  most  commonly  cited 
derivations  of  Lagrange’s  Planetary  Equations  use  Poisson  and  Lagrange  brackets  (see 
Brouwer  and  Clemence),  but  Wiesel  uses  a  simpler  method  using  the  Delaunay  elements 
[27:94],  Avoiding  the  temptation  to  re-derive  the  equations,  the  final  disturbing  function 
form  is  presented  as  [27:98]: 


da  2  dR 
dt  na  dM0 

de  _  1  -  e2  dR  Vl  -  e2  dR 
dt  na2e  dM0  na2e  dco 

di  cot  i  dR  1  dR 

dt  na2Vl  —  e2  dco  na2Vl  —  e2  sin  i  dil 


(100) 

(101) 

(102) 


d.n 

dt 


dR 


na2Vl  —  eA  sin  1 
do)  Vl  —  e2  dR 


di 

cot  i  dR 


dt 


na2e  de  na2Vl  -  e2  di 
dMn  2  dR  1  —  e2  dR 


(103) 


(104) 


(105) 


dt  na  da  na2e  de 
The  disturbing  function  is  given  by  R  in  the  planetary  equations.  Since  the  gravity 
gradient  is  the  largest  source  of  perturbations  in  the  altitude  band  of  concern  in  this  study, 
we  can  limit  the  disturbance  function  to  the  geopotential.  The  zonal  harmonic  given  by 
m  =  0  and  n  =  2,  also  known  as  the  J2  term,  has  the  single  largest  contribution  to  the 
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geopotential  by  at  least  three  orders  of  magnitude  [27:137].  Therefore,  the  disturbing 
function,  R2,  is  approximated  by  the  J2  term  alone.  See  Wiesel  [27:137]  for  its  full 
glorious  form.  After  eliminating  periodic  terms  from  R2,  its  secular  form,  containing  only 
a ,  e,  and  i  is: 


=  ~  2aHl-e^AlSin2(i)  ~  X)  (106) 

Three  of  Lagrange’s  planetary  equations  contain  only  secular  terms:  (103),  (104),  and 
(105).  Upon  substituting  the  partial  derivatives  of  R2sec  in  each  of  these,  the  resulting 
expressions  show  the  approximate  contributions  from  the  J2  disturbing  function  on  the 
fundamental  basis  frequencies  of  a  satellite  orbiting  Earth  (given  in  the  ECI  frame) 
[27:141]: 


co 


.  3Vm/2^®  fh  .  2  \ 
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(107) 

(108) 

(109) 


2.2  Analytical  Mechanics 

The  previous  discussion  of  Keplarian  motion  relied  entirely  on  Newtonian  mechanics 
in  which  the  motion  of  individual  bodies  was  governed  by  vector  quantities  such  as 
momentum  and  external  forces.  Analytical  mechanics  offers  an  alternative  approach  to 
formulate  equations  of  motion  by  way  of  the  whole  system’s  kinetic  energy  and  potential 
energy  [30:45].  Since  these  two  fundamental  quantities  are  scalars,  component-by- 
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component  analysis  with  vector  quantities  is  not  required  in  analytical  mechanics.  This 


departure  from  geometrical  or  physical  coordinates  frees  us  from  a  particular  coordinate 
system  so  that  generalized  coordinates,  q,  can  be  used.  If  the  system  is  characterized  by  n 
generalized  coordinates  and  k  constraints,  the  system’s  degrees  of  freedom  is  given  by  n 
-*[31:31]. 

Two  branches  of  analytical  mechanics  discussed  in  the  proceeding  sections  are 
Lagrangian  Dynamics  and  Hamiltonian  Dynamics.  The  latter  will  be  used  to  generate  the 
equations  of  motion  for  the  topic  of  this  paper,  but  it  cannot  be  done  without  first  wading 
into  Lagrangian  Dynamics  from  which  it  is  formed.  The  primary  distinction  between  the 
two  branches  is  that  Lagrangian  Dynamics  provides  n  second-order  differential  equations 
containing  generalized  coordinates  qfi  =  1,2, ...  ,n)  and  their  velocities  qfi  = 
1,2 ,...,ri),  while  Hamiltonian  Dynamics  provides  2 n  first-order  differential  equations 
containing  generalized  coordinates  qfi  =  1,2,  ...,n)  and  their  momenta  pfi  =  1,2,  ...,n) 
[30:172], 


2.2.1  Lagrangian  Dynamics 

In  1788,  Italian  mathematician  Joseph  Louis  Lagrange  published  his  seminal  work, 
Mechanique  Analytique,  which  contained  the  generalized  Lagrange  Equations  of  Motion 
for  a  mechanical  body  [32:226].  His  famous  equations  were  largely  influenced  by  his 
work  with  Swiss  mathematician  Leonhard  Euler  to  find  an  analytical  solution  to  the 
tautochrone  problem.  Their  use  of  calculus  of  variations  to  find  the  stationary  value  of  a 
definite  integral  resulted  in  the  Euler-Lagrange  differential  equation  [30:57],  Not  by 
coincidence,  the  Euler-Lagrange  equation  is  similar  in  form  to  Lagrange’s  equations  of 
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motion  since  Lagrange  applied  the  same  method  of  calculus  of  variations  and  the 
principle  of  least  action  to  dynamic  systems  to  his  own  equations  [33:587]. 

The  current  form  of  Lagrange’s  equations,  described  here,  was  formulated  by 
William  Rowan  Hamilton  in  1834  and  1835,  when  he  applied  the  principle  of  least  action 
to  a  Lagrangian  function,  L  =  T  -  V,  which  will  be  discussed  below  [34].  To  arrive  at 
Lagrange’s  equations  using  Hamilton’s  method,  we  start  by  defining  a  variation  in  the 
configuration  space  of  generalized  coordinates  to  understand  the  principle  of  virtual  work 
on  a  static  system.  Then  we’ll  use  d’Alembert’s  principle  to  relate  the  principle  of  virtual 
work  to  dynamic  systems.  Finally,  we’ll  use  Hamilton’s  principle,  which  is  simply  the 
integral  of  d’Alembert’s  principle,  to  obtain  Lagrange’s  equations  of  motion. 

2.2.1.1  Variations 

The  basic  problem  of  calculus  of  variations  seeks  to  find  the  path  y(x)  from  ( xi,yi )  to 
(x2,y2)  that  minimize  or  maximize  the  integral  J  =  fX2  f(x,y,  y')  dx  [33:583].  Lagrange 

proposed  small  variations  to  the  dynamical  path  such  that  an  object  can  travel  between 
the  two  points  in  the  same  time  along  any  path.  Consider  the  example  given  in  Figure  8  in 
which  an  object  is  moving  between  points  A(qiA)  and  B(qiB)  in  the  ^-dimensional 
configuration  space.  The  object  at  point  P’  on  a  varied  path  F(q,  q,  t)  is  separated  from 
an  object  at  point  P  on  the  original  dynamical  path  f(q,  q,  t)  by  a  small  variation  Sf  at 
the  exact  same  time,  t.  At  the  endpoints,  the  variation  Sf  is  exactly  zero. 
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Figure  8.  Lagrange’s  variation  in  the  configuration  space. 

Credit:  Vinti  [31:32] 

The  variation  8f  is  called  a  virtual  displacement  and  it  can  vary  any  coordinate  at  a 
fixed  point  in  time  so  long  as  the  variation  is  consistent  with  all  constraints  on  the  system 
[30:53].  Virtual  displacements  are  distinct  from  the  differential  displacement  taking 
place  between  the  time  interval  dt  in  which  forces  and  constraints  may  change  [30:53]. 
Without  delving  deep  into  the  calculus  of  variations  in  our  pursuit  to  develop  Lagrange’s 
equations  of  motion,  we  need  only  show  that  the  operator  8  is  commutable  with  the  d 
designating  differentials  of  displacement.  Start  by  defining  the  virtual  displacement  as 
[31:31]: 

8f  =  F  —  f  (110) 

The  derivative  is  then  [31:31]: 

jtiSf)  =  F-f  (111) 

It  is  also  true  that  [31:31]: 
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(112) 


8f  =  F-f 

Equating  the  variations  from  (111)  and  (1 12)  we  have  [31:31]: 

From  (1 13)  it  is  shown  that  the  variation  of  the  derivative  is  equal  to  the  derivative  of  the 
variation  where /may  be  a  scalar  or  a  derivative.  More  simply,  this  commuting  property 
is  expressed  as  [31:31]: 

Sd  =  d8  (114) 

2.2.1.2  D’Alembert’s  Principle 

Previously,  the  concept  of  virtual  displacement  was  introduced  to  show  the 
commutability  of  the  8  operator,  but  virtual  displacements  can  also  be  used  in  the  same 
way  that  a  differential  displacement  is  used  with  a  force  to  describe  work  on  a  system. 
When  an  applied  force  acts  on  a  system  through  an  infinitesimal  virtual  displacement 
consistent  with  all  system  constraints,  the  work  done  on  the  system  is  zero  [30:60].  This 
is  known  as  the  principle  of  virtual  work  and  is  a  statement  of  static  equilibrium. 
D’Alembert’s  Principle  extends  the  principle  of  virtual  work  to  a  dynamic  system  in 
which  the  time  derivative  of  momenta  are  also  present  [30:65]. 

To  derive  d’Alembert’s  Principle  we  start  with  Newton’s  second  law  for  a  particle  of 
mass,  nik  ,  in  inertial  space  [31:32], 

rnkrk=Fk  +  Ck  (115) 

where  Fk  and  Ck  are  applied  and  constraint  forces,  respectively.  Now  if  we  rearrange  the 
equation  and  apply  the  principle  of  virtual  work  to  the  kth  particle  [30:65]: 
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8W  —  (Ffc  +  Ck  —  mkr fc)  ■  8rk  —  0 


(116) 


The  work  done  by  the  constraint  force  is  zero  since  the  constraint  force  acts  normal  to  the 
virtual  displacement.  So  for  a  system  of  N  particles,  d’Alembert’s  Principle  is  expressed 
mathematically  as  [31:32]: 

N 

m  =  ~  mkrk )  ■  8rk  =  0  (117) 

k= 1 

The  difference  between  the  applied  forces  and  inertial  forces,  Fk  —  mkrk,  is  called  the 
effective  force.  As  such,  d’Alembert’s  Principle  states  that  the  virtual  work  of  effective 
forces  acting  on  a  system  through  an  infinitesimal  virtual  displacement  is  zero. 

For  a  monogenic  system,  d’Alembert’s  Principle  can  be  altered  to  show  the 
relationship  between  generalized  forces  and  generalized  potential  [31:32], 

N 

YFk-Srk  =  -SV(q,t)  (118) 

k= 1 

where  V  is  the  potential.  The  interest  of  this  author  is  its  application  to  satellites  in  which 
it  will  later  be  used  to  represent  the  satellite’s  gravitational  potential  energy. 

2.2.1.3  Hamilton’s  Principle 

Hamilton’s  Principle  is  merely  an  integrated  form  of  d’Alembert’s  Principle,  but  its 
real  distinction  is  the  ease  with  which  Lagrange’s  equations  of  motion  can  be  found.  The 
derivation  of  Lagrange’s  equations  can  be  done  using  either  principle,  but  it  can  be  a 
chore  using  d’Alembert’s  Principle  since  (117)  uses  position  coordinates  that  may  not  all 
be  independent  [30:66],  Hamilton’s  Principle  reduces  the  problem  to  a  scalar  definite 
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integral  so  that  generalized  coordinates  may  be  used,  and  in  the  case  of  a  holonomic 
system,  the  coordinates  are  independent  [30:72], 

Start  by  integrating  (1 17)  from  0  to  t  [31:33]: 


N  t  N  t 

If  F ■  Sv k  dt  —  if  mkrk  ■  8rk  dt 


(119) 


The  integral  on  the  right  can  be  rendered  as  [31:33]: 


fii  _  ff  _  drk  ^  _ 

i  rk  ■  8rkdt  =  Srk  ■ ——dt  =  rk  ■  Srk\  —  rk-dSrk 

o  Jo  dt  0  Jo 


(120) 


From  the  discussion  of  variations  in  2.2. 1.1,  8rk  =  0  at  the  endpoints  of  the  dynamical 
path,  so  the  first  term  on  the  right  goes  to  zero.  Also  recalling  the  commutability  property 

in  which  d(8rk )  =  S(drk)  =  8(rk)dt,  (120)  can  be  rewritten  in  the  form  [31:33]: 


fb  ■  Sfi,  dt  =  -  I  fb  -  Sri,  dt  =  -  -  Sri  dt 


(121) 


Since  the  total  kinetic  energy  for  a  system  of  particles  is  given  by  T  =  ~Y!k=i  mk  the 
right  side  of  (119)  is  just  [31:33]: 


N 

If 


mkrk  ■  drk  dt  =  —  I  ST  dt 


(122) 


For  a  monogenic  system,  (118)  can  be  substituted  into  the  left  side  of  (1 19)  to  get: 


N  i 

If 


Fk  ■  8rk  dt  =  —  I  8V  dt 


(123) 


Upon  substituting  (122)  and  (123)  into  (119),  Hamilton’s  principle  can  be  expressed  as 


[31:33]: 
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8(T  -V)dt  =  0 


(124) 


The  Lagrangian  function  is  defined  as  [31:34]: 


L  =  T(q,q,t )  -  V(q,t ) 


(125) 


So  Hamilton’s  Principle,  in  its  most  generic  form,  is  simply  [31:34]: 


8Ldt  =  0 


(126) 


2.2.1.4  Lagrange’s  Equations  of  Motion 

With  Hamilton’s  Principle  in  hand,  we  can  directly  deduce  Lagrange’s  equations  of 
motion.  Start  by  forming  [31:34]: 


/  dL  dL 


(127) 


The  partial  with  respect  to  t  is  not  included  since  corresponding  points  along  the  varied 
and  dynamical  paths  are  reached  at  the  same  time.  Substituting  (127)  into  (126)  for  the 
case  of  no  constraints: 


V1  ( dL  dL  \ 

Ll{wfqk+Wfqk)dt  =  0 


(128) 


The  second  component  in  parenthesis  can  be  resolved  as  [31:34]: 


ft  dL  Cl  dL  d  dL  f  Cl  d  (dL\ 

{ Wfqkdt  =  L  W*Tt(Sq«)dt  =  wfqk0-l  (129) 


From  the  discussion  of  variations  in  2. 2. 1.1,  8qk  =  0  at  the  endpoints,  so  the  first  term  on 


the  right  goes  to  zero.  Plugging  back  into  (128): 
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(ill©- 

u  k= 1 


dh 


dqk 


8qkdt  —  0 


(130) 


It  follows  from  the  fundamental  lemma  of  calculus  of  variations,  the  coefficient  of  8qk  in 
the  integrand  vanishes  over  the  range  of  0  to  t  [35:34],  It  is  this  coefficient  from  which 
we  deduce  Lagrange ’s  equations  of  motion  [31:34]: 


d  /  dL\ 
dt \dqk)  dqk 


dh 


(131) 


2.2.2  Hamiltonian  Dynamics 

The  equations  of  motion  given  by  Lagrangian  dynamics  consists  of  n  second-order 
differential  equations;  however,  it  is  sometimes  more  convenient  to  express  the  dynamics 
in  terms  of  2 n  first-order  differential  equations  called  Hamilton’s  equations.  Both 
Lagrangian  and  Hamiltonian  dynamics  require  the  use  of  generalized  coordinates,  but 
they  differ  in  the  choice  of  auxiliary  coordinates.  The  Hamiltonian  form  uses  generalized 
momentum,  pk,  rather  than  generalized  velocities,  qk,  as  was  the  case  in  Lagrange’s 
equations.  The  momenta  are  simply  [31:37]: 

dL(q,  q,  t) 


Vk  = 


dqk 


(132) 


Hamilton’s  equations  of  motion  are  given  by  the  time  derivatives  of  the  generalized 
coordinates  and  momenta.  To  derive  them  as  a  function  of  the  variable  set  (q,  p,  t)  rather 
than  ( q,q,t ),  a  Legendre  transformation  must  be  used.  Its  derivation,  shown  in  [30:93], 
leads  to  the  Hamiltonian  function  [31:37]: 
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(133) 


N 


H (q,  p,  t)  =  ^  pkqk  -  L(q,  q,  t) 

k= 1 


The  time  derivative  of  q  can  be  found  by  taking  the  partial  of  H  with  respect  to  the 
momenta  and  recognizing  that  the  q' s  are  a  function  of  q' s  and  p’s  [31:37]: 

N  N 

dH(q,p,t )  .  V  dilk  V  dL  dqk 


dPj 


,  V  adk  sr  OL 

=  qi+Lv^-L^i 


k=1  dVj  £r[dVkdpj 


.  ,  V  /  \ 

=  ^+> 


(134) 


k=l 


=  ij 

For  a  system  with  rectangular  coordinates,  the  summation  term  in  the  second  line  of 

(134)  vanishes  since  the  partial  of  the  Lagrangian  given  by  L  =  Yjk=i^rnk(xl  +Tfc  + 

zjfc)  -  V{x,y,z,t)  is  just  the  momenta  pXk  =  mkxk,  pyk  =  mkyk,  and  pZfc  =  mkzk 

[31:37],  Thus,  Hamilton’s  first  n  equations  of  motion  are  [31:37]: 

dqk  =  dH(q,p,  t)  5 

dt  dpk 

The  time  derivative  of  p  can  be  found  by  taking  the  partial  of  H  with  respect  to  the 
generalized  coordinates  [31:38]: 
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dHjcpp,  t) 
dqj 


N 


-I 

k= 1 


dl 

9qj 


N 


I 

k= 1 


dL  dqk 
dqk  dqj 


-!(*- 


dL_  \  dqk 
dqj  dqj 


dL 

dqj 


dL(<b  t) 

dqj 


(136) 


Substituting  (131)  into  (136)  and  recognizing  the  definition  of  momenta  from  (132) 
[31:38]: 


dH(q,p,t )  d  /  dL\  d  ,  . 

dqj  dt  \dqj  dt 

Thus,  Hamilton’s  second  «  equations  of  motion  are  [31:38]: 

dpk  _  dHjcpp,  t) 
dt  dqk 


(137) 


(138) 


2.3  KAM  Theory 

In  his  famous  1954  address  to  the  International  Congress  of  Mathematicians,  Andrey 
Kolmogorov  first  posed  the  theory  that  a  lightly  perturbed,  conservative,  dynamical 
system  will  exhibit  lasting  quasi-periodic  motion  on  an  invariant  N-torus  [36]. 
Kolmogorov’s  student,  Vladimir  Arnold,  and  German- American  mathematician  Jurgen 
Moser  rigorously  proved  the  theory  for  Hamiltonian  systems  [37;  38].  This  new  approach 
to  stability  problems  in  celestial  mechanics  would  become  the  foundation  for  more  than  a 
half-century  of  work  in  a  field  that  now  bears  their  names:  Kolmogorov- Amold-Moser 
theorem. 
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In  its  simplest  mathematical  representation,  KAM  theory  is  not  unlike  classical 
perturbation  theory  with  the  initial  state  governed  by  [39:6]: 

HE0,<p)  =  h(0  +  £f0,<p)  (139) 

where  Hs  is  a  perturbed  Hamiltonian  which  is  27r-periodic,  h  and  /  are  real-analytic 
functions  representing  an  unperturbed  Hamiltonian  and  a  perturbing  function, 
respectively,  £  is  a  small  («  1)  real  valued  perturbing  parameter,  and  (7,  q>)  are 
symplectic  action-angle  variables  on  the  torus.  One  distinct  difference  between  the  two 
theories  and  proof  of  KAM  theorem’s  value  is  the  speed  at  which  the  solution  is 
converged  upon.  Using  a  sequence  of  canonical  transformations,  the  solution  from 
classical  theory  is  converged  upon  linearly,  if  at  all.  For  example,  in  the  first  step,  the 
initial  Hamiltonian  Hx  =  hx  +  efx  can  be  transformed  to  H2  =  h2  +  £2/2  in  which  the 
order  of  the  perturbation  grows  linearly  to  £2.  On  the  jth  iteration,  the  Hamiltonian  is  of 
the  form  Hj  =  hj  +  £; fj  and  the  perturbation  has  grown  to  a]  [40:43].  With  every 
iteration,  the  denominator  of  the  generating  function  can  grow  arbitrarily  small,  causing  a 
divergence  from  the  solution  with  higher  orders  of  £  [40:39].  This  is  the  well-known 
small  divisor  problem  that  plagued  mathematicians  such  as  Henri  Poincare  [41], 
Kolmogorov’s  theorem  overcame  the  small  divisor  problem  by  converging  upon  the 
solution  quadratically  such  that  after  the  jth  iteration  the  Hamiltonian  is  of  the  form 

Hj  =  hj  +  £  fj.  This  approach  controls  the  small  divisor  in  the  sequence  of  canonical 
transformations  so  that  infinitely  many  iterations  may  be  used  [40:43]. 

This  super-convergent  analytical  technique  was  used  to  prove  KAM  theory  and  show 
that  solutions  to  a  non-degenerate  Hamiltonian  will  persist  on  an  invariant  torus  as  long 
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as  the  perturbations  remain  small — an  entirely  new  approach  to  perturbation  theory. 
Classical  perturbation  theory  seeks  to  approximate  the  solutions  and  then  explore  its 
evolution/stability  from  a  fixed  initial  condition,  whereas  KAM  theory  doesn’t  concern 
itself  with  the  motion  incurred  from  preassigned  initial  conditions,  but  instead  explores 
the  dynamic  stability  in  phase  space  using  a  set  of  fixed  frequencies  that  govern  quasi- 
periodic  motion  [40:40],  In  the  integrable  case  when  £  =  0,  the  phase  space  solution  will 
lie  on  an  invariant  torus  with  a  set  of  N  fundamental  frequencies.  When  £  is  sufficiently 
small  and  the  frequencies  are  sufficiently  incommensurate  (satisfying  the  diophantine 
inequality),  a  solution  is  quickly  converged  upon  that  remains  on  the  invariant  torus  (a 
condition  of  stability).  As  £  grows,  the  torus  is  deformed  or  displaced  until  it  ceases  to 
exist  [40:42], 

2.3.1  Torus  Visualization 

To  fully  understand  the  nature  of  the  torus,  its  dimensionality  must  be  explored.  For 
the  typical  Hamiltonian  system  encountered  in  an  earth  satellite  orbit,  the  simplified 
point-mass  experiences  three-degrees  of  freedom  in  three-dimensional  “native”  space 
characterized  by  the  generalized  coordinates  q  =  (qlt  q2,  q2)  and  their  conjugate 
momenta  p  =  (p i,P2,Ps)-  Later,  we  will  develop  the  exact  Hamiltonian,  H(q,  p), 
represented  by  these  variables.  In  the  meantime,  it  is  enough  to  say  that  KAM  theory 
seeks  to  map  this  lightly  perturbed  Hamiltonian  to  one  represented  by  new  coordinates, 
Q  =  (Qi.Qz'Q-i)’  and  momenta,  P  =  (P1(  P2,  P3),  in  which  only  the  momenta  appear  in 
the  new  Hamiltonian,  76  (P).  (See  Wiesel  [19:7]  for  further  discussion  of  the  new 
Hamiltonian  and  generating  function  as  approximated  with  the  Delaunay  variables.)  The 
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absence  of  the  new  coordinates  merely  implies  from  Hamilton-Jacobi  theory  that  the 
coordinates  have  constant  frequencies.  The  three  present  momenta  combined  with  the 
three  missing  coordinates  produce  6-dimensions  defining  the  phase  space  of  the  torus. 
This  is  equivalent  to  three  sets  of  action-angle  pairs,  each  pair  defining  a  circuit  on  the 
torus.  The  new  coordinates  are  the  angles  which,  like  their  native  counterparts,  enable 
three-degrees  of  freedom.  The  three  new  momenta  are  integrals  of  motion  which  mandate 
that  the  solutions  lie  on  a  three-dimensional  manifold  which  is  topologically  equivalent  to 
a  three-torus  [42], 

Any  torus  larger  than  a  two-torus  is  impossible  to  intuitively  comprehend  since  an  n- 
torus  must  exist  in  at  least  n  +  1  dimensions.  The  three-torus  of  interest  to  this  work  is 
impossible  to  visualize  since  it  is  given  by  six  dimensions,  but  no  loss  is  incurred  since  it 
can  still  be  communicated  perfectly  in  mathematical  form.  It  is  possible  to  consider 
partial  visualization  techniques  with  lower  order  tori.  For  example,  a  one-torus  (given  by 
one  action-angle  pair)  with  one  degree  of  freedom  exhibits  bounded  motion  that  looks 
like  a  circle  (two-dimensional).  A  two-torus  (given  by  two  action-angle  pairs)  with  two 
degrees  of  freedom  exhibits  bounded  motion  with  respect  to  two  closed  circuits  that 
resemble  a  doughnut  (three-dimensional).  For  higher  order  tori,  an  approach  similar  to 
Poincare  mapping  can  be  used  to  visualize  the  orbit  in  a  lower  dimensional  subspace 
[27:25].  In  this  case,  a  two-torus  can  be  generated  by  taking  the  “cross-section”  of  the 
three-torus  which  effectively  divorces  one  of  the  coordinates  and  its  conjugate  momenta 
from  the  others.  Figure  9  depicts  a  standard  two-torus. 
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Figure  9.  Cross-sectional  view  of  an  invariant  two-torus  defined  with  the  typical  action 
angle  variables.  In  this  case  the  new  momenta  resemble  the  actions  which  define  the 
shape  of  the  torus  onto  which  the  motion  may  be  projected  and  the  angles  give  the 

position. 
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III.  Methodology 


Orbital  data  from  the  ISS  will  be  fit  to  a  phase  space  trajectory  that  exhibits  drifting 
toroidal  motion  over  a  dense  continuum  of  adjacent  tori  due  to  stochastic  and 
deterministic  forces.  This  is  distinct  from  multiply  periodic  motion  on  a  single  torus  that 
accounts  for  purely  deterministic  forces.  The  fitting  process  begins  with  the  construction 
of  an  initial  reference  torus  from  a  numerically  integrated  orbit  that  includes  the  full 
geopotential  up  to  order  and  degree  m,n  =  20.  Wiesel  has  shown  that  deterministic  tori, 
such  as  the  reference  torus,  exhibit  exquisite  accuracy  to  within  a  few  meters  over  a 
decade  when  compared  to  the  integrated  orbit  (considered  as  the  truth  source).  The  next 
step  is  to  account  for  stochastic  effects  in  the  orbital  data  that  are  not  modeled  by  the 
reference  torus.  Intuitively,  the  stochastic  forces  impart  an  offset  to  the  reference  state 
vector.  This  appears  as  a  displacement  in  the  reference  torus’  initial  momenta  and 
coordinates,  the  impact  of  which  is  motion  on  an  adjacent  torus.  Bayesian  estimation  is 
used  to  find  a  vector  offset  that  minimizes  the  residual  between  the  actual  orbit  and  the 
reference  torus.  Sequential  outputs  from  the  Bayes  filter  provide  updates  to  the  reference 
torus  that  should  produce  a  very  good  estimate  of  the  ISS  orbit  over  an  extended  period 
of  time.  The  entire  process  of  constructing  the  reference  torus  and  fitting  ISS  data  to  the 
reference  torus  is  explored  and  developed  here. 

3.1  Reference  KAM  Torus 

Unfortunately,  there  doesn’t  exist  an  easy  homeomorphism  from  the  three- 
dimensional  space  (given  by  generalized  coordinates,  qt,  and  their  momenta,  pt)  to  the 
six-dimensional  phase  space  of  a  torus  (given  by  actions,  /*,  and  angles,  <p*  ). 
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Nevertheless,  the  common  periodic  nature  of  the  physical  coordinates  and  the  torus  angle 
variables  makes  an  N-tuple  Fourier  series  a  possible  solution  for  mapping  between  both 
spaces.  There  are  generally  two  approaches  used  for  constructing  the  torus  Fourier  series: 
1)  trajectory-following  techniques  that  perform  a  Fourier  decomposition  of  data  from 
long  numerically  integrated  orbits,  and  2)  iterative  techniques  that  find  successively 
better  approximations  of  the  series  using  Hamilton’s  equations  [43:147].  Binney  and 
Spergel  pioneered  the  first  approach  for  non-integrable  galactic  dynamics  in  1982  and  it 
has  since  become  the  most  straightforward  and  robust  approach  for  orbits  [43:149]. 
Wiesel  has  recently  shown  two  variants  of  this  type  for  earth  orbits  [17;  18].  The  first  is  a 
one-pass  approach  that  uses  least  squares  to  fit  a  Fourier  series  to  the  integrated  orbit  and 
the  second  is  a  two-pass  approach  that  identifies  the  fundamental  frequencies  on  the  first 
pass  and  then  seeks  the  series  coefficients  on  the  second  pass.  The  latter  will  be  employed 
later  in  this  section. 

The  aim  of  both  methods  is  to  construct  the  following  finite  Fourier  series  truncated 
to  order  M  =  Qw  hltmit,  hlimit), 

M 

^  Dj(l)exp(ij-(p )  (140) 

j=-M 

where  Dj  are  the  complex  series  coefficients,  j  is  the  index  summation  vector,  and  j  ■  (p 
are  the  associated  frequency  combinations.  The  more  conventional  real  form  gives  the 
Fourier  series  as: 

M 

9(0  =  C(0jo . o)N  +  ^  ft  cos(/  '  Q)  +  sin(/  ■  Q)}  (141> 

j  =—M 
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(142) 


Q(t)  =  a(t  -  to)  +  Q(t0) 

Here  q(t)  are  the  physical  trajectory  coordinates  in  the  time  domain,  Q(t)  are  the  torus 
angle  variables  incremented  linearly  in  time,  fl  is  a  set  of  N  fundamental  or  basis 
frequencies,  Cj  and  Sj  are  the  Fourier  coefficients  for  each  coordinate,  and  j  is  the  index 
summation  vector.  The  basis  frequencies  describe  the  satellite’s  underlying  periodic 
behavior  and  are  easily  extracted  from  a  power  spectral  density  (PSD)  plot  of  the  phase- 
space  trajectory.  The  Fourier  coefficients  are  just  the  amplitudes  of  the  Fourier  transform 
at  each  integer  combination  of  the  basis  frequencies.  The  index  summation  is  an 
incrementing  scheme  that  allows  for  integer  combinations  of  the  basis  frequencies.  It  also 
ensures  that  the  first  non-zero  term  in  j  isn’t  negative  so  as  to  avoid  repeated  angle 
combinations  in  the  numerical  routine,  the  effects  of  which  are  not  lost  to  symmetrical 
trigonometric  properties,  specifically  sin(0)  =  —  sin(0)  and  cos(0)  =  cos(— 0). 

A  variation  of  Jacques  Laskar’s  Numerical  Algorithm  of  the  Fundamental  Frequency 
(NAFF)  will  be  the  machinery  of  choice  for  finding  frequencies  and  coefficients  to 
construct  the  torus  Fourier  series.  The  NAFF  is  a  technique  that  approximates  the 
truncated,  continuous  Fourier  transform  (TCFT)  so  that  prominent  spectral  lines  can  be 
identified  without  the  destructive  effects  of  aliasing  and  leakage  which  will  be  attended  to 
momentarily  [44;  45],  The  frequency  approximations  at  these  peaks  can  then  be  used  in  a 
least  squares  approach  to  identify  the  basis  set.  With  the  basis  set,  Laskar’s  NAFF  would 
typically  be  used  to  extract  the  coefficients  from  individual  spectral  peaks;  however,  the 
spectral  decomposition  method  used  here  is  a  variation  of  the  NAFF  that  extracts 
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coefficients  from  clusters  of  peaks  using  the  analytical  form  of  the  truncated,  continuous, 
Fourier  transform  (ATCFT)  in  the  fitting  process  [46]. 

One  is  left  to  ponder  the  exact  relationship  between  the  Fourier  series  parameters  and 
the  action-angle  variables  that  describe  the  satellite  position  on  the  torus.  From  classical 
mechanics,  the  constant  actions  define  the  shape  of  invariant  tori,  while  the  angles  are  the 
coordinates  on  the  tori.  Mathematically,  the  canonical  angle  coordinates  are  those  given 
in  (142),  but  their  conjugate  action  momenta  are  nowhere  to  be  found  (explicitly)  in  the 
Fourier  series.  The  momenta  represent  the  dimensions  in  phase  space  that  are  directed 
away  from  the  torus  surface  and  are  only  implicitly  present  through  the  series 
coefficients.  Nevertheless,  they  can  still  be  calculated  explicitly  from  the  Poincare 
integral  invariants  [17;  18;  19]: 


where  Qj  are  the  angle  coordinates,  p*  and  qt  are  the  native  coordinates  and  momenta, 
and  r \  is  a  fundamental  contour  about  the  torus.  As  the  system  oscillates  around  the  torus, 
the  time  derivative  of  the  coordinates  is  equivalent  to  the  basis  frequencies.  Since  the 
action  momenta  are  constant  on  a  torus,  their  time  derivative  is  simply  zero,  which  infers 
the  system’s  Hamiltonian  function  is  only  a  function  of  action  momenta.  This  makes 
sense  because  the  Hamiltonian  is  conserved.  As  such,  the  Hamiltonian  equations  of 
motion  are  resolved  as: 


dQi  _  dJC(P) 
dt  dPi 


Cli(P) 


(144) 
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(145) 


dPi  _  dK(P ) 

dt  dQi 

From  (144)  it  is  apparent  that  the  frequencies  in  the  orbital  motion  are  a  function  of  the 
action  momenta.  Thus,  the  influence  of  the  action  momenta  on  the  Fourier  series  is 
apparent  in  the  Fourier  coefficients  as  amplitudes  of  the  frequencies  and  their 
combinations  in  time. 

The  subsections  that  follow  will  document  the  exact  procedures  for  generating  the 
reference  torus,  but  a  preliminary  roadmap  is  given  here: 

•  Define  the  satellite’s  dynamics  with  a  Hamiltonian. 

•  Numerically  integrate  the  Hamiltonian  using  the  satellite’s  reference  position 
and  velocity  as  the  initial  condition. 

•  Convert  the  integrated  trajectory  from  the  time  domain  to  the  frequency 
domain  using  a  finite  Fourier  transform. 

•  Extract  the  basis  frequencies  and  Fourier  coefficients  using  Laskar’s  method. 

•  Construct  a  Fourier  series  representation  of  the  KAM  torus. 

3.1.1  Earth  Satellite  Dynamics 

A  fixed  torus  about  the  Earth  can  only  be  constructed  from  a  lightly  perturbed 
autonomous  Hamiltonian  system.  The  reference  satellite  orbit  includes  perturbations 
from  the  geopotential  and  ignores  all  nonconservative  forces,  so  the  only  temporal 
variations  are  induced  from  the  rotating  Earth.  To  freeze  the  geopotential  and  obtain  an 
autonomous  dynamical  system,  the  ECEF  frame  is  imposed  on  the  system.  The 
Hamiltonian  dynamics  can  then  be  formulated  as  done  by  Wiesel  [17:1-2], 

Start  by  defining  the  generalized  coordinates  in  the  ECEF  frame  and  their  inertial 
velocities  resolved  in  the  ECEF  frame: 
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(146) 

The  Lagrangian  function  given  by  (125)  is  no  longer  a  function  of  time  since  the 
dynamics  are  autonomous  in  the  non-rotating  frame, 

L  =  T(q,q)-V(q) 


(147) 


=  -  ((x  -  cu0y)2  +  (y  +  w0x)2  +  z2)  —  V (, q ) 


where  T (q)  is  the  expanded  geopotential  from  (95).  The  canonical  momenta  are 
calculated  from  (132)  and  are  identically  equal  to  the  inertial  velocity  components, 


V  = 


x  -  (x )©y 
y  +  co0x 
z 


(148) 


where  o>0  is  the  Earth’s  rotation  rate  in  canonical  units  of  0.0588335998  rad/TU  where  1 
time  unit  (TU)  is  13.446852  minutes.  Substituting  the  appropriate  terms  into  the 
Hamiltonian  function  given  by  (133)  yields: 

1 

H(q,  p,t)=-  (pi  +Py  +  pi)  +  M@(yVx  ~  XPy ) 


oo  n 

P  V  V  (  r 


LL  [rz  p”m(cos« 


(149) 


n= 0  77i=0 


X  \Cnm  COS  (ml)  +  Snm  sin(mA)] 

/?0  is  the  equatorial  radius  of  the  Earth  equal  to  6378.137  km  or,  in  canonical  units,  1 
distance  unit  (DU);  p  is  the  Earth’s  gravitational  constant  equal  to  398600.4418  km3/s2 
or,  in  canonical  units,  1  DU3/TU2;  and  as  before,  n  and  m  are  the  degree  and  order  of  the 
geopotential  expansion,  P™  are  the  Legendre  polynomials,  and  Cnm  and  Snm  are 
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dimensionless  coefficients  from  the  gravity  model.  Finally,  the  radius,  r,  the  geocentric 


latitude,  8,  and  the  east  longitude,  A,  are  found  from: 

r  =  yjx2  +y2  +  z2 

7 

sin  8  = 


V*2  + 


r 


y 

tan  A  =  — 
x 


(150) 

(151) 

(152) 


The  Hamiltonian  equations  of  motion  are  determined  by  taking  the  partials  of  the 
Hamiltonian  function  as  specified  in  (135)  and  (138)  such  that. 


Px  3" 

Py  —  co^qx 


Q 

P- 


Pz 

dV(q) 

a,®Py  dqx 
dV(q) 

aJ®Px  dqy 
dV 


(153) 


L  dqz  J 

where  the  partials  of  V(q)  are  the  geopotential  force  upon  the  satellite  and  are  only 
dependent  on  the  satellite’s  position.  Given  the  autonomous  nature  of  the  Hamiltonian 
function,  the  equations  of  motion  are  independent  of  time  and  the  Hamiltonian  is  a 
constant  of  motion. 


3.1.2  Numerical  Integration 

A  Hamming  fourth-order  predictor-corrector  algorithm  is  used  to  integrate  the 
satellite’s  Hamiltonian  equations  of  motion  across  a  one-year  timespan  from  time  -T  to 
time  T  to  produce  coordinate  data,  /(t).  This  amounts  to  a  6-months  forward  and  6- 
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months  backward  integration  with  initial  conditions  given  at  time  0.  The  deliberately 
symmetric  forward  and  backward  integration  is  convenient  for  a  finite  Fourier  transform, 
but  more  importantly,  it  limits  the  accumulation  of  total  error  compared  to  a  1-year 
forward  integration  to  time  2T  or  a  one  year  backward  integration  to  time  -2T.  Since  the 
associated  local  error  induced  by  the  algorithm  is  of  order  h5  where  h  is  the  time-step,  the 
steps  are  limited  to  40.340556  seconds  for  a  total  of  400,000  steps  in  each  direction. 

The  total  error  growth  can  be  conveniently  checked  at  each  time-step  by  calculating 
the  absolute  change  in  the  Hamiltonian  function: 

AH  =  \H(t)  —  H(t0)\  (154) 

The  Hamiltonian  should  be  a  constant  of  motion,  yet  A H  =£  0  because  the  integrator  is 
not  symplectic.  In  other  words,  the  time  evolution  of  the  equations  of  motion  do  not 
possess  a  conserved  Hamiltonian.  This  is  a  common  result  for  most  numerical  integration 
schemes  since  the  integration  does  not  conserve  the  symplectic  two-form,  dp  A  dq. 

The  numerical  integration  of  the  ISS  trajectory  is  initialized  at, 

Tol  /— 0.682606200715413X  /  0.339437086117151  vf 
_  =  -0.082702006542665  ,  -0.892456547165198  (155) 

LpJ  l\  0.797396182426368  /  V  0.197152638998895  /. 

where  the  position  vector  is  given  in  DU  and  the  momenta  vector  is  given  in  DU/TU.  All 
geopotential  terms  up  to  order  and  degree  m,n  =  20  are  included  in  the  integration.  Figure 
10  depicts  the  first  30  days  of  the  orbit  propagation  in  the  TDR  frame.  By  comparison,  a 
full  one-year  representation  of  the  orbit  would  cover  virtually  the  entire  globe  between  52 
degrees  north  and  52  degrees  south  of  the  equator. 
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Figure  10.  30-day  representation  of  ISS  orbit  propagation 


As  seen  from  Figure  11,  the  maximum  amplitude  error  is  roughly  10-12  which  is 
only  four  orders  of  magnitude  from  double  precision  accuracy.  Given  that  the 
Hamiltonian  and  state  vector  are  scaled  to  order  unity,  this  miniscule  error  assures  that 
the  propagated  Hamiltonian  is  only  faintly  perturbed  from  the  original  one  and  KAM 
theorem  will  still  apply. 
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Figure  11.  Hamiltonian  error  from  Hamming  fourth-order  integrator. 

3.1.3  Fourier  Transform  and  Spectral  Analysis 

The  transformation  of  the  numerically  integrated  coordinate  data,  /(t),  from  the  time 
domain  to  the  frequency  domain  is  done  using  a  finite  Fourier  transform  over  the 
symmetric  time  interval  [-T,  T\.  The  typical  Fourier  transform  of  a  function  over  an 
infinite  time  span  is  done  by, 

n  OO 

T(v)  =  I  /  (t)e~2nvit  dt  (156) 

J  —  OO 

where  v  is  the  cycle  frequency.  The  Fourier  transform  assumes  infinite  periodicity  in  the 
signal,  but  over  an  arbitrary,  finite  time  interval  it  is  doubtful  that  the  signal  endpoints  are 
of  the  same  value.  The  consequence  of  this  discontinuity  is  a  phenomenon  known  as 
spectral  leakage  in  the  Fourier  transform  [47],  To  inhibit  this,  the  signal  is  multiplied  by 
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a  window  function  that  forces  the  signal  to  start  and  end  at  zero  amplitude.  Laskar’s 


NAFF  uses  a  Hanning  window  of  the  form  [44,  45], 


t\\p 


(157) 


where  p  is  the  order  of  the  cosine  function  and  T  is  the  frequency  interval. 

The  window  function  is  shown  graphically  for  various  values  of  p  in  Figure  12  for 
timespan  [-7,  7].  With  the  window  function,  the  finite  Fourier  transform  over  the 
timespan  [ -7 ,  7]  becomes: 

't' 


T(v) 


(158) 

-  M /M*-2*-*  if)  *  (f)  dt 

The  domain  of  the  Fourier  transform  can  be  expressed  as  cycle  frequency,  v,  or  angular 
frequency,  oj  =  2nv,  but  the  latter  will  typically  be  used  for  the  torus. 
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Figure  12.  Hanning  window  of  order  p. 


Choosing  an  appropriate  window  power  is  essential  to  dissecting  spectral  content 
from  the  transformed  data.  The  effect  of  increasing  the  window  power  is  shown  in  Figure 
13  by  applying  a  Fourier  transform  on  a  single  spectral  line  of  unit  amplitude  over  the 
timespan  [-1,  1]  as  demonstrated  by  Laskar  [44:10]  and  Wiesel  [18:4]: 


(159) 


In  this  example,  the  signal  is  a  gate  function  and  has  a  single  spectral  line  at  o>0  =  0,  but 
Figure  13  shows  sidelobe  oscillations  of  a)T  =  n/T  which  appear  from  the  cosine  term  in 
(157).  As  p  increases,  the  main  lobe  of  the  a>0  spectral  line  is  broadened  and  the  sidelobes 
around  it  fall  off  more  rapidly.  The  amplitude  of  the  spectral  line  remains  the  same  with 
increasing  p,  so  the  advantage  of  higher  order  window  functions  is  an  accelerated 
convergence  upon  the  basis  frequency.  One  disadvantage  is  found  in  signals  containing 
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Figure  13.  Graph  of  sidelobe  effects  from  Tq(oj)  to  T4(oj). 

integer  combinations  of  frequency  sets  in  which  one  frequency  is  significantly  smaller 
than  the  other: 

3  (ftm,ftn)  e  ft  I  ftm  «  ftn,m  *  n  (160) 

where  mathematical  jargon  compresses  the  conditional  statement  to  read:  if  there  exists 
(3)  a  pair  of  frequencies  (ftm,  ftn)  as  an  element  of  (6)  the  basis  set  ft,  such  that  ( I )  ftm 
is  significantly  smaller  than  ftn  where  m  and  n  are  not  equal.  When  this  occurs,  such  as 
the  case  of  cascading  harmonics  from  the  small  apsidal  regression  frequency,  the  higher 
order  window  functions  can  “swallow”  nearby  spectral  content.  This  is  known  as  spectral 
shadowing. 

Ultimately,  the  choice  of  the  Hanning  window  power  is  contingent  upon  the  signal 
characteristics,  but  Laskar  shows  for  KAM  solutions  that  the  accuracy  of  the  frequency 
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analysis  (order  1  IT  2p+1)  generally  improves  by  increasing  p.  He  suggests  using  the 
highest  possible  value  of  p  until  precision  begins  to  decrease.  Wiesel  has  typically  used  a 
Hanning  window  of  order  p= 2  to  avoid  spectral  shadowing.  This  was  found  to  be 
sufficient  for  the  integrated  ISS  data. 

Once  the  data  is  represented  in  the  frequency  domain  with  an  appropriate  window 
function,  the  spectrum  must  be  decomposed  into  a  set  of  precise  basis  frequencies  and 
Fourier  coefficients.  Wiesel  does  this  by  specifying  an  approximate  set  of  basis 
frequencies  (described  in  §3. 1.3.1)  followed  by  a  Newton-Rhapson  algorithm  that  seeks 
the  maximum  spectral  power,  T  =  \(p\2,  in  the  vicinity  of  integer  combinations  of  the 
approximate  basis  frequencies  [17:5].  As  the  peaks  are  identified,  the  basis  set  can  be 
determined  from  a  least  squares  solution  of  the  frequency  approximations.  With  the 
newly  estimated  basis  set,  the  coefficients  of  the  Fourier  transform  are  determined 
according  to  the  following: 

6(o,o,  ...,o)w  =  9^(0)  (161) 

Cj  =  23Fh(o)j)  (162) 

Sj  =  -234>(w;)  (163) 

is  the  complex  Fourier  transform  at  the  composite  frequency  <x>j  =  /■  ft,  or  in  the 
case  of  the  constant  term,  C(0  0  0yv,  at  co  =  0. 

3.1.3.1  Basis  Frequency  Approximations 

The  starting  guess  for  the  three  basis  frequencies  in  the  TDR  frame  are  determined 
from  motion  in  the  two-body  orbit  with  the  addition  of  secular  rates  from  the  J2 
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disturbing  function  (see  the  Lagrange  Planetary  Equations  from  §2. 1.6.1).  The  largest 
frequency  is  the  anomalistic  frequency  found  in  the  orbit  plane.  In  terms  of  the  classical 
orbital  elements,  it  is  approximately  the  Keplarian  frequency  with  the  mean  anomaly 
correction: 


n  ~  H  3Vm/2#0  (l  .  2  r  \ 

1  ~  Va3  2a7/2(l  -  e2)3/2  V2  Sm  W  1 ) 


(164) 


The  second  largest  frequency  appears  in  the  true  of  date  equator.  It  is  approximately  the 
Earth’s  rotation  rate  (which  must  be  negative  since  the  node  appears  to  move  clockwise 
about  the  TDR  z-axis)  with  the  nodal  regression  correction: 


Ll2  ~  tn©  —  _  „  , - ^ttcos  i 


(165) 


2a7/2(l  —  e2)2 

The  smallest  frequency  is  approximately  the  apsidal  regression  rate  that  characterizes  the 
rotation  of  the  line  of  apsides  due  to  the  Earth’s  oblateness: 


flo  ~  — 


2a7/2(l  —  e2)2 


(^sin2(0  —  2) 


(166) 


3.2  Motion  Near  an  Earth  Satellite  KAM  Torus 

The  motion  of  an  Earth  orbiting  satellite  is  dependent  on  more  than  just  conservative 
gravitational  forces,  so  the  reference  torus  alone  will  not  suffice  for  most  satellites 
(especially  those  in  low  Earth  orbit).  Wiesel  claims  that  the  satellite’s  physical  state 
vector  can  be  modeled  with  a  contribution  to  the  reference  torus  from  non-deterministic 
forces  such  that, 

X(t)  =  [x,  y,  z,  vx,  vy,  vz]  =  Xref(t)  +  SXst0C(t)  (167) 
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where  8Xst0C(t)  is  the  non-deterministic  or  stochastic  displacement.  The  stochastic 
displacement  can  occur  parallel  to  the  torus  surface  and  perpendicular  to  it.  The  parallel 
portion  results  in  a  displacement  to  the  torus  coordinates,  Q,  which  can  be  accounted  for 
on  the  surface  of  the  reference  torus.  The  perpendicular  portion  cannot  be  accounted  for 
on  the  surface  of  the  reference  torus  since  any  displacement  to  the  torus  momenta,  P, 
steps  off  the  torus  surface.  Thus  the  stochastic  displacement  must  be  broken  into  two 
parts:  that  within  the  reference  torus  and  that  beyond  the  reference  torus.  Reflecting  both 
parts,  the  displacement  can  easily  be  resolved  in  terms  of  the  torus  coordinates  and 
momenta  using  the  chain  rule  [48:2], 


-  -  ,  N  dX(t)  dY(t)  -  N 

X(t)  =  XreAt)  +  ^  SY(t0) 

f  dY  (t)  av(to) 

y  m  ,  dm  dm  am  dm 

-  Xre/(t)  +  ~^—^8Y{t0)  +  x  N  oY(t0) 


dQ(t)dY(t0) 


dP(t)  dY(t0) 


(168) 
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Parallel  Perpendicular 

where  the  local  displacement  with  respect  to  the  reference  torus  is  given  by: 


SY(t)  = 


SQ(t ) 
8P(t) 


(169) 


The  process  for  determining  the  local  displacements  is  left  for  the  next  section,  but 
first,  let  us  direct  our  attention  to  how  they  are  used  for  both  parallel  and  perpendicular 
motion.  There  are  two  ways  to  account  for  the  parallel  motion  labeled  in  (168).  Not 
surprisingly,  the  most  direct  method  simply  adds  the  parallel  8X  component  to  the  torus 
generated  Xref,  exactly  as  indicated.  The  preferred,  indirect  method  adds  5Q  to  the  torus 
coordinates  which  makes  it  unnecessary  to  add  the  parallel  8X  component  to  Xref.  Either 
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way,  the  contributions  of  5Q  must  be  accumulated  over  time  and  stored  in  Q(t0)  to 
minimize  error  growth  in  the  reference  torus  and  avoid  a  breakdown  of  the  linearization 
model.  Unlike  parallel  motion,  there  is  only  one  way  to  account  for  the  perpendicular  SX 
component;  it  must  be  added  directly  to  Xref  since  P  remains  unchanged  on  the  reference 
torus. 

Even  though  the  momenta  do  not  change  on  the  reference  torus,  the  actual  orbit  is 
characterized  by  new  momenta  formed  by  the  displacement,  P  +  8P,  which  is  still 
constant,  but  describing  an  adjacent  torus  with  different  frequencies.  This  frequency  shift 
must  be  accounted  for  in  the  reference  torus  coordinates  if  it  is  to  accurately  approximate 
motion  nearby.  Recall  that  the  state  vector  for  the  torus  coordinates  and  momenta  is  given 
by: 


Y(t)  = 


<2(0 

ft  ■  (t  -  to)  +  Q(t0) 

P(fX 

constant 

(170) 


It  was  previously  shown  in  (144)  that  the  constant  frequencies,  ft,  are  a  function  of  the 
momenta,  so  if  the  momenta  change  by  an  incremental  amount,  so  too  do  the  frequencies 
[48:7]: 


5Q  =  5ft  =  ^5P  (171) 

dP 

Since  the  torus  coordinates,  Q,  increment  linearly  with  time  as  seen  in  (170),  a  slightly 
different  frequency  manifests  as  a  linear  drift  in  the  corresponding  torus  coordinate.  We 
directly  observe  the  drift  magnitude  by  integrating  (171)  to  get: 
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When  the  offset,  5Q,  is  added  to  the  reference  torus  coordinates,  it  should  minimize  the 
error  between  the  actual  orbit  (potentially  on  an  adjacent  torus)  and  the  reference  torus 
(likely  nearby).  That  is  to  say,  the  8 Q  update  to  the  reference  torus  should  maintain  a 
close  proximity  to  the  actual  orbit  by  accounting  for  parallel  local  motion  near  the 
reference  torus.  With  the  update,  the  reference  torus  coordinates  are  simply: 

Q(t)=a-(t-to)  +  Q(to)  +  8Q  (1V3) 

Wiesel  indicates  that  the  reference  torus  initial  coordinates  must  be  routinely  updated 
with  the  offset,  5Q,  to  remain  within  the  region  for  which  the  linear  observation  model  is 
a  valid  approximation  of  the  satellite’s  motion. 

Next,  let  us  direct  our  attention  to  the  partial  derivative  matrices  from  (168)  which 

will  be  needed  for  sequential-batch  estimation.  5Q(t)/dY(t0)  is  shown  to  be  the 
following  3x6  matrix: 


3Q(t)/dY(t0)  = 


8Q  dQ 


dQo  9P0 


(174) 


Taking  the  partials  of  (173)  with  respect  to  Q0,  it  is  clear  that  the  left  Jacobian  is  a  3x3 
identity  matrix.  The  right  Jacobian  is  much  more  complicated.  If  (144)  is  expanded  with 
the  chain  rule  to  include  a  small  change  in  the  fundamental  frequencies  [49], 


dQ  -  dil  ^  -  d2X  ^ 

-P  =  ft  +  —  5P  =  Sl  +  ^^8P 
dt  dP  dP2 


(175) 


5Q/dP0  can  then  be  approximated  from  the  second  partial  derivative  of  the  Hamiltonian 


function,  7C[19:8]: 


X  = 


r 


2P2 


-  P2(0<£  + 


fh RUPj  -  3 pf) 


(176) 


4  P13P35 

where  the  first  term  is  from  two-body  motion,  the  second  is  a  vestige  of  the  partial 
derivative  of  the  generating  function  with  respect  to  time  (an  artifact  of  using  the  TDR 
frame  for  the  torus  coordinates),  and  the  last  term  is  the  J2  potential  in  terms  of  the 
Delaunay  momenta.  The  first  partial  of  X  is  the  frequency  given  by  (144),  or  written  as  a 
scalar  [48:7]: 

dX 


dPi 


(177) 


The  second  partial  derivative  of  X  is  then  the  typical  Jacobian  [48:7]: 


d2X 
dP2  ~ 

Applying  the  finishing  touch  to  (174),  the  full 


10  0 


dQ(t)/dY(t0)  = 


0  1 


0 


0  0  1 


d£l 

dP 


matrix  is 

d01 

At 

dPi 

dnx 

At 

dP2 

301  ■ 

At 
dP3 

df 12 
dPx  At 

d02 

dfl2 

At 

dfl3 
dPj.  At 

d03 

dfl3 
a/’s At 

where  the  matrix  components  are  [48:7]: 

ggi  _  3/r2  3/r4/2f?2(p|-p2) 

dPt  P*  P15P35 

oa±  =  on2  =  9h%R®p2 

dP2  dP1  2P*Pf 

ggi  djh  9/t4;2P|(P32  -  5 Pj) 

dP3  dP1  P*P* 


(178) 


(179) 


(180) 

(181) 

(182) 
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(183) 


an2  3jj%Rl 

dP2  2P13P35 

an2  =  an3  =  ls^hR^ 

dP3  dP2  2P13P36 

an3  3^/2fl|  (2P32  -  15Pf) 

dP3  2  P13P37 

In  a  similar  manner,  dP(t)/dY(t0)  is  easily  shown  to  be: 


(184) 

(185) 


dP(t)/dY(t0)  = 


dP  dP 


(186) 


L5q0  ap0J 

From  (144),  the  momenta  are  constant,  so  the  left  Jacobians  is  a  3x3  zero  matrix  and  the 
right  Jacobian  is  a  3x3  identity  matrix. 


dP(t)/3Y(t0)  = 


rO  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


(187) 


The  first  three  rows  of  the  6x3  Jacobian  matrix,  dX(t) / dQ (t),  are  all  that  is  required 
to  produce  displacements  to  the  native  coordinates  and  can  be  determined  directly  from 
the  partial  derivatives  of  the  torus  Fourier  series  previously  specified  by  (141)  such  that, 

M 

dq(t)/dQ;(t)  =  ^  ji{-Cj  sin(/-  Q)  +  5y  cos (/■  Q)}  (188) 

j=-M 

A  less  accurate  alternative  is  to  approximate  the  full  Jacobian  matrix  from  the  derivatives 
of  the  linearized  two-body  problem.  This  is  the  same  method  required  for  determining 


aX(t)/dP(t). 


It  follows  that  the  6x6  observation  partial  derivative  matrix,  dX(t)/dY(t),  is  just  the 
concatenation  of  5X(t)/dQ(t)  and  5X(t)/dP(t).  The  first  three  columns  of 
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5X(t)/dY(t)  are  5X(t)/dQ(t)  and  the  last  three  columns  are  5X(t)/dP(t).  The  matrix, 
3X(t) / dY (t),  cannot  be  determined  directly,  but  it  can  be  formulated  from  a  mapping 

between  two  separate  Jacobians,  namely  dY/ dZ  and  d Z/ dX,  where  Z  is  the  vector 
formulation  of  the  classical  orbital  elements: 


Z(t)  =  [M,D.,(o,a,e,i]T 


(189) 


This  mapping  is  possible  since  the  canonical  torus  coordinates  are  very  nearly  equivalent 
to  the  first  three  classical  orbital  elements.  As  a  consequence,  we  can  write  Z  as  [48:2]: 

Z(t)  =  [Q1,Q2,Q3,a,e,i]T  (190) 


The  first  partial  derivative,  dY/dZ,  is  the  Delaunay  elements  Jacobian.  Wiesel  has 
shown  that  upon  ignoring  Earth’s  rotation  rate  and  limiting  the  gravity  field  to  a  single 
point  mass  term,  as  in  the  two-body  problem,  the  torus  coordinates  and  momenta  revert  to 
the  Delaunay  elements  given  in  (84)  -  (89)  [19:7].  Thus,  the  state  vector  Y  can  be  written 


in  terms  of  the  Delaunay  elements  as  a  reasonable  approximation  [48:6]: 


Y(t)  = 


Qi  =  M 
Q2  =  0, 

Q3  =  0) 

Pi  = 

P2  =  yfjias/ 1  —  e2  cos  i 
_  P3  =  1  —  e2  . 


(191) 


The  Delaunay  elements  Jacobian  is  discerned  in  the  usual  fashion: 
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3Y/dZ  = 


1  0 
0  1 
0  0 

0  0 


0  0 


0  0 


0 

0 

o  - 

0 

0 

0 

0 

0 

0 

a  px 

dP1 

d  Px 

da 

de 

di 

dP2 

dP2 

d?2 

da 

de 

di 

dPs 

dPs 

dPs 

da 

de 

di  - 

where  the  matrix  components  are  [48:6]: 


d?2 

da 

dPz 

de 

d_?2_ 

di 


dP1  1  fju 
da  2\a 


dP1 

de 

dP1 

di 


=  0 


=  0 


i  m 


—  —v  1  —  e2  cos  i 
2^a 


jua 

—e  \- - tCos  i 

1  —  e2 


—yfjlay/ 1  —  e2  sin  i 


(192) 


(193) 

(194) 

(195) 

(196) 

(197) 

(198) 

(199) 

(200) 

(201) 
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The  second  partial  derivative,  dZ/ dX,  is  the  two-body  elements  Jacobian  in  which 
the  classical  orbital  elements  are  linearized  in  the  TDR  frame.  Start  by  characterizing  the 
Jacobian  matrix  as: 


dZ/dX  = 


~dM 

dM 

dM 

dx 

dy 

dz 

dCL 

an 

an 

dx 

dy 

dz 

dot 

dco 

dco 

dx 

dy 

dz 

da 

da 

da 

dx 

dy 

dz 

de 

de 

de 

dx 

dy 

dz 

di 

di 

di 

dx 

dy 

dz 

dM 

dM 

dM 

dvx 

dvy 

dvz 

an 

an 

an 

dvx 

dVy 

dvz 

dco 

da) 

da) 

dvx 

dvy 

dvz 

da 

da 

da 

dvx 

dvy 

dvz 

de 

de 

de 

dvx 

dvy 

dvz 

di 

di 

di 

dvx 

dVy 

dvz 

(202) 


Given  the  definition  of  the  semimajor  axis  in  (41),  the  partial  derivative  components  in 


the  fourth  row  of  (202)  are: 


da  ii4x 
dx  8 a2r3 

da  jU4y 
dy  8 a2r3 

da  /t4z 
dz  8 a2r3 

da  [i3vx 
dvx  8a2 

da  n3vy 
dvy  8a2 


(203) 

(204) 

(205) 

(206) 

(207) 
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(208) 


da  /j.3vz 


dvz  8a2 

Given  the  definition  of  the  magnitude  of  the  eccentricity  vector  in  (69),  the  partial 
derivative  components  in  the  fifth  row  of  (202)  are  [48:3]: 


de  If  dex  dey  de ZN 

6x~dff+ey~dff  +  ez 


dx  ||e|| 
1 


dx 


de 

dy 

de 

dz 

de 


de 


X 


de. 


|e||  yGx  dy  +  £y  dy 


+  e 7 


1 

W\ 


dex  dey 

ex  ~~z  I-  ey  ~ ~  h  e 

dz  y  dz 


de. 


dvx 

He|| 1 

[xdvx 

de 

1  | 

(  dex 

dVy 

l|e|| ' 

V  dVy 

de 

1  j 

[e  — 

dvz 

l|e|| 1 

{  Xdvz 

+  e. 


de 


dv 


de z\ 
dy) 

dez\ 
z  dz  ) 

y-  +  ez- ^ 


dv 


Xj 


+  e. 


dey  dez 

- H  ez - 


dv,, 


dv 


yy 


dey  dez 

4-  p - u  p  - 

+  dvz  +  dVzy 


(209) 

(210) 

(211) 

(212) 

(213) 

(214) 


The  partial  derivatives  of  the  eccentricity  vector  required  in  (209)  -  (214)  are  determined 
from  (68)  [48:3]: 


dex 

=  UV2. 

ll  \IX 

dx 

dex 

dey 

1  /fixy 

dy 

dx 

jU  v  r3 

dex 

_  dez  _ 

1  mxz 

dz 

dx 

rr3 

dey 

dy 


M  My2 

r  o 

y i  J 


VXVy^ 

vxvz ) 


(215) 

(216) 

(217) 

(218) 
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dey  dez  1  myz  \ 

dz  ~  dy  ~  /jlr3  VyVz) 

(219) 

dez  1/  n  \iz2  > 

oz  fjL\  r  r6  j 

)  (220) 

dVx= /  ^  z,o 

(221) 

Qj  Qj 

VT  * 

II 

“fc  1 

* 

1 

£ 

(222) 

dex  _  1 

(2xi?z  zvx) 

dvz  \i 

(223) 

II 

1 

>? 

1 

X 

£ 

(224) 

II 

"fc  1 

T 

X 

£ 

1 

N 

(225) 

Qj  Qj 

N  ^ 

II 

1  m 

'To' 

Vi 

1 

N 

(226) 

dez  1 

=  (2zvx  xvz) 
dvx  n 

(227) 

(228) 

dez  1  ,  \ 

dvz  =  ^XV‘  yv^> 

(229) 

Given  the  quadrant  conclusive  definition  of  inclination  in  (71),  the  partial  derivative 
components  in  the  sixth  row  of  (202)  are  [48:4], 


di 

dx 


-1 


'1  dhz  hz 


dhy 


dhx 


2  \  h  dx  h3  [hx  dx  +  dx 


+  h7 


dh7 


dx  \ 


(230) 
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di 

dy 

di 

dz 

di 

dvx 

di 

dVy 

di 

dv7 


— 1  /ia/iz  hz 


dhx  dhv  dh7 

hy  — h  hV  — - h  Hn 


,u  \2  \h  dy  h3  rx  dy  '  "y  dy  '  "z  dy  J 

!-(f) 


— 1  fldhz  hz 


-(I) 


dhx  dhv  dhz 
h  — -  +  h  — -  +  h  — - 
tu  dz  h3  L  az  37  az  z  azj 


■l  /ia/iz 


fh7\2\hdvx  b3 

!-(f) 


—  1  fldhz  hz 


dhx  dhy 

hxd^  +  hyd^  +  hz 


dh7 


dvY 


-m 


h7\2\hdvy  h3 


dhx  dhv  dhz 

avy  ^  avy  ayy 


— i  /i  a/iz  /iz 


-®) 


^  n2  \hdv7  h 3 


dh7 


dhY  dhy 

hx~d77  +  hy~d77  +  hzdv7 


(231) 


(232) 


(233) 


(234) 


(235) 


where  the  specific  angular  momentum  vector,  h,  is  given  in  (13)  and  its  partial 
derivatives  are  [48:4]: 


ah/ax  = 


dhx 

dhx 

dhx 

dhx 

dhx 

dh 

-1 

dx 

dy 

dz 

dvx 

dVy 

dvz 

dhy 

dhy 

dhy 

dhy 

dhy 

dh 

y 

dx 

dy 

dz 

dvx 

dVy 

dvz 

dhz 

dhz 

dhz 

dhz 

dhz 

dhz 

dx 

dy 

dz 

dvx 

dVy 

dvz_ 

"  0 

vz 

1 2y 

0 

—z 

y 

-vz 

0 

vx 

z 

0 

—x 

Vy 

-vx 

0 

-y 

X 

0 

Given  the  quadrant  conclusive  definition  of  RAAN  in  (75),  the  partial  derivative 
components  in  the  second  row  of  (202)  are  [48:4]: 


(236) 
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da 


ky 


dh. 


dx  h%  +  hj,  \  x  dx 


da 


ky 


dh. 


dy  h2  +  h2\  dy 


da 


dz  h%  + 


da 


ky 


dhy 

dz 


dh. 


—  ky 


-h 


-h 


dhx 

dx 

dh x 
y  dy 

dhx 
y  dz 

dhY 


-  -  h  — --h 

dvv  hi  +  hi  \  x  dvY  y  dvY 


da 


x  ,Lxx,Ly 

1 


ky 


dhy 


~  hy, 


dhy 


dVy  h%  +  h%  \  x  dvy  y  dvy/ 


da 


hy 


dhy 


dvz  h2  +  h2\  x  dvz 


dhx 

dv7 


(237) 

(238) 

(239) 

(240) 

(241) 

(242) 


Given  the  definition  of  argument  of  perigee  in  (76),  the  partial  derivative  components  in 
the  third  row  of  (202)  are  [48:5]: 


doi 

dx 


-1 


x  _  (Key  -  hM2  ^  +  ^ 
V  ejh2  +  h2  ) 


'  dhx  dev  dhv 
e  — -+h  — —  —  e  — - 

ydx+  x  dx  x  dx 


hy  dx)  (hx6y  hy6^ 


i—3 


dey 


V  hx  +  hy  V  dx 


dey  dez' 
+  6y  ~fa+ez  ~dl, 


(243) 


+ 


e{h2x  +  h2) 


(  dhx  dhy' 

^\hx~d^  +  hy~d^. 
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doo 

dy 


-1 


dh-y  de , 


dhx 


x  (cVSFT7iJ  v  y  ay  ay 
V  eJhx  +  hy  ) 


de 


hy  Qy  )  (jlxey  h yex ) 


f  dex  dey  de z> 

ey  —  h  ~  I-  e7 


y/h%  +  hj\  x  dy  y  dy  z  dy 


1  (  dhx  dhy 

- \ - 777 1  hx  — 1-  hy  — — 

e{hi  +  hjf2  \  dy  gy 


doi 

dz 


-1 


dhy 


de v  dhy, 


1  - 


^Z^{eM^yVy  +ft‘  »*  a* 

eJhZTh*  ) 


hy  ^  (h-x^y  hyex ) 


j-3 


3ex  dey 


3e7 


[V^TT^V  ^ 


+  ey  dz  +  6z  az 


+ 


e{hl  +  hj) 


dhx  dh 
hx  -r — I-  h 


3/2 1'^  az 


5z 


a  00 

a^v 


-l 


dhx  dey  dhy 

ey  — - h  /ir  - e y 


/,  ,  \  2  le,//i2  +  /i2  V  y  avx  x  avx  xavx 

1  Mxey-/iyex\  (eV^x^y  \  xxx 

V  ejhl  +  /i|  j 


ae 


hy  gv  )  {hxey  hye^ 


(  dex  dey  dez 

6v  ~  I-  £y  t  H  — 

x  Hn  017 


V^x  +  hj\  x  dvx  ydvx 


e(/i|  +  hj) 


(  dhx  dhy 

^r*a7  +  /lya7 


(244) 


(245) 


(246) 
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doo 

dvv 


-1 


i  _  fhxey-hyex\2  leV^  +  Vy 

V  ejwthi  ) 


dhx  dev  dhv 

e  — -  +  h  — -  — e  — - 

y  dv  x  dv  x  dv 
y  \  u  uy  u  vy  u  vy 


dv  )  (.hxey  hyex) 


y/h$  +  hj\  xdv. 


dex  dey  de7\  (247) 

+  ev  — - 1-  e7 


X  ^y  y  dvy 


dvv 


e{hx  +  hj ) 


/  dhx  dhy N 

^[hxd^r  +  hyd^:. 


Vy 


d  0) 
dv7 


-1 


dhy 


1  _  ( hxey  -  hyex\2 1 eV^  +  h* 


Gy  dv7  +  hx  dv. 


dev  dhy, 


z  x  dvz 


ey/W+h^J 
~  ^  dv  )  _  (hxey  ~  hyex) 


,-3 


dey 


[yfWThjK  ^VZ 


dey  dez 
+  e  — —  +  e  — - 

y  dv7  z  dv7 , 


(248) 


1  (  dhx  dhY\ 

+  e{h*  +  h^f2  \hx~d^z+hy~d^) 

where  the  components  of  the  angular  momentum  vector,  eccentricity  vector,  and  their 
partials  are  as  before.  Since  the  argument  of  perigee  requires  a  quadrant  check,  the 
partials  must  be  multiplied  by  -1  if  the  “z”  component  of  the  eccentricity  vector  is  less 
than  zero,  ez  <  0. 


The  final  set  of  partial  derivatives  in  the  first  row  of  (202)  is  the  most  complicated. 


Starting  with  the  definition  of  the  mean  anomaly  at  epoch  (53)  [48:6], 


dM 

dx 


,  dE  de 

(1  —  e  cos  E)— - —sin  E 

dx  dx 


(249) 


dM  dE  de 

—  =  (1  -  ecosE)- - —sin E  (250) 

dy  dy  dy 
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(251) 


dM  dE  de 

-r—  =  (1  -  e  cos  E)  - —  sin  E 

dz  dz  dz 


dM 

dvx 


dE  de 

(1  —  e  cos  E)  — - - — sin  E 

dvx  dvx 


(252) 


dM 

dvy 


dE  de 

(1  —  e  cos  E)  — - - —  sin  E 

dVy  dVy 


(253) 


dM  dE  de 

- — =  (1  -  e  cos  E)  - - — sin^  (254) 

dvz  dvz  dvz 


Recall  that  the  partial  derivatives  of  the  eccentricity  are  given  by  (209)  -  (214);  however, 


the  partial  derivatives  of  the  eccentric  anomaly  come  from  (66)  such  that  [48:5]: 


sec2 


Qdv 

die 


1  1  —  e  \  de 

1  +  e  (1  +  e)2J  dx 


sec2 


v\  dv 
2 )  dy 


1  1  —  e  \  de 

1  +  e  (1  +  e)2J  dy 


sec2 


Qdv 
dz 


/  1  1  —  e  \  de 

\1  +  e  (1  +  e)2/  dz 


(255) 


(256) 


(257) 
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dE 

dvx 


dE 

dVy 


i+rrl^n2g)[2^ 


1  —  e 


1  +  e 


■sec 


Qdv 
dvY 


1  / v\  1  +  e  (  1  1  —  e  \de 

"I*3"©, 


1+irltan2©i2'J 


1  —  e 


+  e  (1  +  e)2Jdvx 


Qdv 


- sec 

l  +  e  Wdv* 
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~  2 tan  \2.) ,  +  (TTeV) 


+  e  (1  +  e)2/  dvy 


(258) 


(259) 


dE 

dvz 


1  + 


1  —  e 
l  +  e 


tan2  g) 


1  —  e 


l  +  e 


sec 


2(V)+L 

^2/  dvz 


l  +  e 


1  —  e 


d  +  e 


1  —  e  \  de 
(1  +  e)2J  dvz 


(260) 


The  partial  derivatives  of  the  eccentric  anomaly  introduce  the  need  for  the  partial 
derivatives  of  the  true  anomaly,  v.  From  (81)  we  know  that  for  the  non-degenerate  case. 


V  =  0)  —  u 


dv  dco  du 
dX  ~  dX  dX 


(261) 

(262) 


Since  the  partial  derivatives  of  argument  of  perigee  have  already  been  determined,  all  that 
is  required  are  the  partial  derivatives  of  the  argument  of  latitude,  u,  given  by  (80)  [48:5]: 
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With  the  completion  of  (202),  the  observation  partial  derivative  matrix,  dX/dY,  can 
be  formed  by  inverting  the  product  of  the  Delaunay  elements  Jacobian  and  the  two-body 
elements  Jacobian  [48:6]: 


dX 

dY 


( dY  dZ\ 
\dZ  dX.) 


(269) 


One  method  to  verify  the  two-body  analytical  formation  of  dX/ dY  is  to  numerically 
differentiate  a  collocation  formula.  Here  I’ve  chosen  Stirling’s  formula  to  get  [50:112], 

( SYn0 )  =  ^  (xm_2  -  8xm_1  +  8xmi  -  xm2 )  (270) 

where  xm  are  native  coordinates  and  momenta  determined  from  the  torus  Fourier  series  at 
one  observation  time.  For  each  observation,  (270)  must  be  evaluated  with  five  data 
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points,  {xrn-2>xrn-i’Xm0>xm.1>xm2)m_1  6-  Each  is  formed  from  a  slight  variation  to  the 
core  8Yno  such  that  ( 8Yn_2 , 8Yn  i,  8Yno,  8YUl,  8Yn2)n=1  are  all  separated  by  a  small  step 
size,  h.  In  other  words,  (8YnQ  —  2 h,8YnQ  —  h,8YnQ,8YnQ  +  h,8Yno  +  2h)  .  During 

the  least  squares  fitting  process,  which  will  be  discussed  in  the  next  section,  h  can  be  set 
to  any  small  step  size  on  the  first  iteration  since  8YUq  =  0.  On  subsequent  iterations  and 
for  Bayesian  filtration,  it  is  preferred  that  the  step  size  be  scaled  relative  to  the  magnitude 


of  8Yno.  The  overall  6x6  partial  derivative  matrix  is  then  formed  according  to  the  typical 


Jacobian, 


dX/dYi  = 
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(271) 


3.3  Toroidal  Sequential-Batch  Estimation 

The  process  of  estimating  an  actual  orbit  with  a  reference  torus  is  a  new  one.  With 
both  slow  and  fast  variational  forces  acting  on  a  satellite,  uncertainties  in  the  dynamics, 
and  errors  in  observational  data,  a  stochastic  approach  must  be  applied  to  determine  a 
best  estimate  of  the  state.  The  primary  question  arising  from  this  concerns  when  and  how 
state  updates,  5Y0,  should  be  incorporated  into  the  reference  torus.  Should  it  be  done  at 
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each  successive  observation  time  (as  done  with  Kalman  filters)  or  should  it  be  done  at  the 
beginning  of  a  batch  of  observations  over  an  extended  period  of  time  (as  done  with 
nonlinear  least  squares)?  If  the  answer  is  the  latter,  how  should  the  state  and  covariance 
be  brought  forward  in  time  for  every  new  batch  of  data?  Without  a  previously  published 
or  accepted  technique  for  filtering  the  data  and  estimating  the  orbit  with  a  KAM  torus, 
this  author  considered  all  possibilities  from  the  traditional  carte  du  jour  for  nonlinear 
systems. 

In  the  Kalman  filter  class  of  differential  correction  techniques,  two  were  considered: 
linearized  and  extended  Kalman  filters.  Given  the  complexities  of  generating  a  reference 
torus  at  this  time,  the  so-called  extended  Kalman  filter  was  immediately  discounted  since 
it  requires  re-integrating  the  orbit  (thus  generating  a  new  torus)  at  each  observation  time 
[24:736].  The  linearized  Kalman  filter  which  avoids  having  to  re-integrate  the  reference 
trajectory  also  will  not  work  since  it  does  not  utilize  state  updates  in  subsequent 
predictions  [24:735].  It  is  crucially  important  that  routine  updates  be  provided  to  the 
reference  torus  to  prevent  the  actual  orbit  from  drifting  beyond  a  region  correctable  by 
linearized  models.  The  usual  sequential-batch  least  squares  method,  sometimes  called  a 
Bayes  Filter,  was  also  dismissed  for  its  inability  to  propagate  the  state  and  covariance 
through  time  [24:716].  Nevertheless,  a  variation  of  the  Bayes  filter  that  does  provide 
updates  to  the  reference  torus  was  deemed  appropriate. 

Our  version  of  the  Bayes  filter  produces  an  update  to  the  initial  phase  angles,  Q0,  of 
the  torus  after  every  batch  of  observational  data  is  processed.  Subsequent  torus 
predictions  are  produced  from  the  updated  epoch  coordinates.  In  addition  to  phase  angle 
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updates,  Bayesian  estimation  demands  a  priori  data  passed  from  batch  to  batch.  On  the 
first  batch,  the  nonlinear  least  squares  algorithm  is  used  to  generate  the  first  phase  angle 
update  and  covariance  matrix.  But  before  either  the  nonlinear  least  squares  algorithm  or 
the  Bayes  filter,  the  data  must  be  “pre-filtered”  to  eliminate  rogue  data.  Thus  we  begin 
with  the  development  of  the  so-called  pre-filter. 


3.3.1  Pre-filtration 

A  pre-filtration  routine  was  developed  by  this  author  to  eliminate  fallacious 
observational  data,  Z,  that  appears  as  improbable  spikes  in  the  residual  vector  given  by: 


G  = 


Xi  — 

yi-9i 


(272) 


lzi-Zi\ 

where  the  predicted  coordinates  from  the  state  vector,  X,  are  notated  by  the  “hat.”  The 
residual  from  GPS  data  is  generally  smooth  except  in  the  case  of  a  maneuver,  but  even 
post-maneuver  the  residual  is  still  quite  smooth,  just  with  a  different  magnitude.  As  such, 
outlier  data  is  easy  to  identify.  The  entirety  of  the  data  is  not  pre-filtered  all  at  once; 
rather,  it  is  pre-filtered  in  batches  in  sequence  with  the  nonlinear  least  squares  and  bayes 
filter.  This  ensures  the  residuals  between  the  observation  data  and  the  reference  torus 
remain  small  compared  to  outlier  data. 

The  pre-filter  evaluates  the  data  on  three  passes  using  three  different  techniques  to 
identify  rogue  data.  The  first  pass  merely  eliminates  data  outside  of  four  standard 
deviations.  The  second  pass  generates  a  moving  average  and  moving  standard  deviation 
from  eleven  data  points  -  five  before  and  five  after  the  current  observation,  except  at  the 
beginning  and  end  of  the  data  file  which  uses  the  first  and  last  eleven  data  points. 
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respectively.  Data  outside  of  3  moving  standard  deviations  from  the  moving  average  are 
eliminated.  The  third  pass  evaluates  the  absolute  value  of  the  slope  between  the  current 
and  previous  data  point  and  if  the  slope  exceeds  10  km/TU,  the  current  data  is  rejected. 
The  first  data  point  must  be  good  for  the  third  pass  to  work.  The  author  monitored  which 
data  were  rejected  and  the  total  effect  of  the  three  passes  was  flawless  for  an  entire  week 
of  ISS  data  sampled  at  one  minute  intervals. 


3.3.2  Nonlinear  Least  Squares  (NLS) 

Once  the  data  is  pre-filtered,  the  iterative  non-linear  least  squares  algorithm  is 
initiated.  At  zero  epoch,  the  residual  between  the  observation  data  and  the  reference  torus 
should  be  small.  If  it  is  not,  the  reference  torus  is  probably  not  sufficient  for  an  accurate 
estimate  and  must  be  refined.  When  the  residual  is  small  enough  (less  than  1  m  is  ideal), 
the  state  vector  can  be  predicted  at  each  observation  time  using  a  null  displacement  on 
the  first  iteration: 


5K0  =  [0  0  0  0  0  Of  (273) 

Next,  determine  the  residuals  from  observation  data,  Z,  and  estimated  data,  X,  using 
(272).  The  observation  data  of  the  ISS  is  from  the  Global  Positioning  System  (GPS) 
which  shows  very  little  noise.  We  approximate  the  total  instrumental  covariance  matrix: 

10 
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(274) 


If  both  position  and  velocity  were  of  interest  to  this  study,  the  observer  would  be  defined 
as  a  6x6  identity  matrix;  however,  since  only  the  native  position  coordinates  are  of 
interest,  define  the  observer  as: 
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1  0  0  0  0  O' 

H=  010000 
.0  0  1  0  0  0. 


(275) 


From  the  local  dynamics  given  in  (168),  form  the  6x6  state  transition  matrix  as: 


dX(t)  dY (t) 
dY (t)  dY (t0) 


ax(t)  dQ(t)  dX(t )  dP(t) 
aQ(t)aY(t0)  3P(t)  dY(to) 


(276) 


where  the  partial  derivatives  are  the  same  as  those  developed  previously.  For  each 


observation,  form  the  3x6  version  of  the  observation  partial  derivatives  matrix  according 


to: 


Ti  =  H[<P(ti,t0)]  (277) 

Since  the  observational  data  is  pre-filtered  to  eliminate  fallacious  points,  rejection 
processing  during  the  nonlinear  least  squares  algorithm  is  turned  off;  nevertheless,  if  the 
reference  torus  predictions  are  accurate  enough,  each  observation  can  be  processed  to 
reject  data  in  which  the  residual  is  beyond  3  standard  deviations  of  the  instrumental 
noise: 


Qx,y,z  =>  rejected  (278) 

If  the  data  is  not  rejected,  add  to  the  running  sums  of  the  inverted  covariance  matrix,  P-1, 
and  the  vector,  p: 


p-i  =  £rr 

i 

(279) 

P  =  ^  !■[  Q-'r, 

(280) 

l 
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For  each  position  coordinate,  add  the  square  of  the  residual  to  a  running  sum.  Once 


all  observations  have  been  processed,  calculate  the  RMS  of  residuals  according  to: 


RMSXiyiZ  = 


N 


'Ltr? 


i  'x,y,z 


N 


(281) 


Check  for  convergence.  If  the  RMS  of  residuals  for  each  coordinate  is  less  than  the 
standard  deviation  of  the  instrument  noise,  the  solution  is  at  hand. 


RMSX:y  Z  < 


converged 


(282) 


If  the  solution  has  not  converged,  improve  the  estimate  of  SV0  at  epoch  by  adding  the 
product  of  the  covariance  matrix  and  the  vector,  p: 

SY0  (+)  =  SY0  (-)  +  Pp  (283) 

Then  repeat  the  process  from  the  beginning  using  the  new  estimate  of  SY0.  Note  that  we 

are  not  calculating  SY0  at  every  observation;  we  are  doing  so  for  the  whole  batch  of  data. 
Also,  this  update  will  be  discarded  if  the  convergence  criteria  are  met  in  the  current 
iteration. 

Once  converged  or  when  the  maximum  number  of  iterations  has  been  reached,  log 
(5T0  and  the  observation  covariance  matrix  to  pass  forward  to  the  Bayes  filter  for  the 
second  batch  of  data. 


3.3.3  Bayesian  Filtration 

Begin  by  pre-filtering  the  data  as  before,  but  predictions  from  the  reference  torus  are 

determined  using  5F0  from  the  least  squares  routine  (only  if  on  the  second  batch)  or  from 
the  previous  Bayes  filtration.  To  avoid  confusion  with  the  new  displacement  that  will  be 
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determined,  we  will  call  this  SY0old.  Next,  bring  the  initial  phase  angles  and  covariance 

matrix  forward  to  the  current  epoch  using  SY0old.  The  update  to  Q0  is  determined  from 
(172)  such  that: 

Qonew  =  Qoold  +  SQ(t0new,  t0old)  (284) 

where  the  value  of  SQ  is  accumulated  from  the  old  epoch  to  the  new  epoch.  The 
covariance  matrix  is  propagated  using  the  state  transition  matrix,  ®(tonew>toold),  from 
(276)  and  using  the  covariance  matrix,  P0id,  from  the  previous  epoch: 

PneW(~)  =  *Poid<*>T  (285) 

With  the  phase  angle  update  and  a  priori  data  now  available  at  the  new  epoch,  begin 
the  first  iteration  using  the  previous  correction  to  the  torus  state  vector: 

«?< w(-)  =  «?oow  (286) 

As  before,  calculate  the  residuals  and  observation  partial  derivative  matrix  (277)  for  each 
observation,  and  accumulate  the  summation  matrices  (279)  and  (280).  After  all 
observations  have  been  processed,  the  inverse  covariance  matrix  is  given  by: 

Pnew(+ )  =  Pnew(~)  +  ^  Q^T j  (287) 

i 

Next,  check  for  convergence  using  the  RMS  for  each  coordinate  as  presented 
previously  in  (281)  and  (282).  If  the  solution  has  not  converged,  improve  the  estimate  of 
<SK0  at  epoch  according  to: 

«?»_(+)  =  PneA+KPnU-^O^i-)  ~  ?0„» J  +  P)  <288) 
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Since  ft  is  constant  and  Q0  is  constant  for  a  single  batch,  S Q  and  5P  are  all  that  change. 


so: 

Y0  (— )  —  Y0  =SY0  (-)-$?„  (289) 

unew v  y  unew  unew v  y  unew  v  y 

As  with  the  nonlinear  least  squares,  repeat  the  process  with  every  new  estimate  of  <ST0 
until  the  solution  converges  or  until  the  maximum  number  of  iterations  has  been  reached. 

Then  log  SY0  and  the  observation  covariance  matrix  to  pass  forward  to  the  next  batch  of 
data. 


3.4  Stochastic  Orbit  Modeling 

Modeling  stochastic  orbits  requires  discrete  knowledge  of  both  gravitational  and  non- 
gravitational  perturbations.  The  reference  torus  allows  us  to  distinguish  between  them 
since  it  already  contains  the  earth’s  geopotential  effects.  The  added  dynamics  near  the 
reference  torus,  as  estimated  by  the  least  squares  fitting  process  in  the  previous  section,  is 
a  result  of  non-gravitational  perturbations  and  can  be  parameterized  for  apparent 
stochastic  behavior.  Stochastic,  by  definition,  infers  a  random  nature,  but  because  it  is 
believed  motion  near  the  torus  exhibits  an  ostensible  orderliness  that  can  be  modeled  for 
a  period  of  time  under  similar  environmental  conditions,  this  author  shall  refer  to  torus 
perturbations  as  pseudo-stochastic. 

To  produce  stochastic  predictions  via  the  reference  torus,  the  initial  phase  angles,  Q0, 
must  be  varied  with  time  (think  of  this  as  parallel  motion  on  the  surface  of  the  torus)  as 

do  the  momenta  offsets,  SP,  which  allow  the  satellite  to  move  off  the  surface  of  the 
reference  torus  (think  of  this  as  perpendicular  motion ).  The  behavior  of  the  initial  phase 
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angles  and  momenta  offsets  are  first  determined  through  Bayesian  estimation  which  may 
be  sensitive  to  the  number  of  observations  included  in  each  batch.  During  the  estimation 
process,  the  torus’  initial  phase  angles  are  updated  after  every  batch  which  encapsulates 
the  influence  of  non-gravitational  perturbations  and  the  momenta  offsets  are  stored  for 
later  use.  Assuming  these  empirical  parameters  are  well-behaved  and  can  be  modeled 
with  a  time-varying  curve,  their  propagation  can  provide  autonomous  updates  to  the 
reference  torus  for  non-gravitational  perturbations.  For  now  on  this  author  will  refer  to 
these  as  the  pseudo-stochastic  parameters.  Given  a  reference  torus  and  the  pseudo¬ 
stochastic  parameters,  one  should  be  able  to  predict  the  non-chaotic  orbit  of  any  satellite. 

3.5  Observational  Data 

The  procedures  outlined  in  this  chapter  may  be  applied  to  any  satellite  in  a  non- 
chaotic  earth  orbit.  As  such,  virtually  any  satellite  could  be  used  to  verify  and  validate  the 
procedures  since  it  is  not  the  practice  of  the  launch  community  to  place  satellites  in 
chaotic  orbits.  The  ISS  was  chosen  as  the  test  case  in  this  study  for  two  reasons.  First,  it 
is  crucial  that  the  observational  data  be  extensive  and  well  documented  to  eliminate  or 
isolate  external  variables  in  the  study.  Since  human  spaceflight  missions  have  the  highest 
priority  for  tracking  and  predictions,  the  ISS  offers  extensive  data  options  as  well  as 
highly  specialized  trajectory  personnel  that  monitor  the  orbit  fulltime.  The  ISS  is  also  an 
ideal  vessel  for  studying  stochastic  predictions  since  it  is  the  largest  manmade  object 
orbiting  the  earth.  In  its  LEO  orbit,  air  drag  imparts  a  significant  perturbing  force  on  the 
vast  surface  areas  spanning  the  ISS.  The  current  configuration  of  the  ISS  measures  357 
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feet,  end-to-end  [51].  Figure  14  shows  an  excellent  juxtaposition  of  the  ISS  with  the 


famous  Lambeau  Field  which  measures  360  feet  including  the  end  zones. 


Figure  14.  The  International  Space  Station  measures  357  ft,  end-to-end,  whereas  an  NFL 
football  field  measures  360  ft.  Credit:  T.  Brian  Jones 


Observational  data  was  provided  by  the  ISS  Trajectory  Operations  and  Planning 
Branch  (DM33),  Mission  Operations  Directorate,  Johnson  Space  Center  in  Houston, 
Texas.  The  data  was  transmitted  from  the  ISS’s  primary  onboard  GPS  receiver/processor. 
The  parameters  used  include: 
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Table  1.  Onboard  Primary  GPS  Parameters. 


LADP06MDQ000H 

x  position,  TDR  coordinate  system  [m) 

LADP06MDQ001H 

y  position,  TDR  coordinate  system  [m) 

LADP06MDQ002H 

z  position,  TDR  coordinate  system  (m) 

LADP06MDQ003R 

x  velocity,  TDR  coordinate  system  [m/s) 

LADP06MDQ004R 

y  velocity,  TDR  coordinate  system  [m/s) 

LADP06MDQ005R 

z  velocity,  TDR  coordinate  system  [m/s) 

LADP06MDQ006W 

Nav  time  [milliseconds  since  the  beginning  of  the 
week  which  starts  on  Sunday  at  midnight) 

LADP06MDQ452W 

GPS  time  [seconds  since  Jan  6, 1980  at  00:00  UTC) 

LADP06MD2946W 

State  quality  [l=good  or  0=poor) 

The  data  was  sampled  at  1 -minute  intervals  from  00:01:00.000  UTC  on  25  March 
2010  to  12:18:01.000  UTC  on  23  April  2010.  This  corresponds  to  42,000  observations. 
To  give  a  sense  of  the  orbit’s  coverage,  the  pre-filtered  data  appears  in  Figure  15. 


Figure  15.  Pre-filtered  ISS  data  from  25  March  -  23  April  2010. 
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IV.  Results 


Recalling  back  to  §1.6,  three  primary  topics  were  to  be  explored  by  this  thesis.  The 
first  two  were  addressed  in  the  previous  chapter  with  a  derivation  of  the  linearized 
equations  of  motion  near  tori  followed  by  their  application  to  an  estimation  routine  for 
fitting  observational  data  to  a  continuum  of  tori.  This  chapter  will  answer  the  final  and 
most  important  question  of  whether  or  not  stochastic  predictions  can  be  generated  from  a 
reference  torus  and,  if  so,  for  how  long.  The  chapter  begins  with  results  from  the 
reference  torus  construction  machinery  followed  by  a  simple  verification  of  the  Bayesian 
filtration  algorithms.  The  chapter  concludes  with  the  estimation  and  prediction  results. 

4.1  Reference  Torus  Construction 

A  relatively  good  fit  was  obtained  between  the  reference  torus  and  the  integrated 
orbit,  the  results  of  which  will  be  discussed  shorty,  but  first,  we  begin  with  the  difficulties 
surrounding  its  development.  The  basis  frequencies  proved  to  be  the  biggest  challenge. 
For  the  ISS,  the  approximate  basis  set  from  Keplarian  and  perturbation  theory,  equations 
(164)  -  (166),  were  found  to  be: 


Table  2.  J2  basis  frequencies  for  spectral  analysis. 


Coordinate 

Frequency  (rad/TU) 

9.24464666422975e-001 

11 2 

-5.96726955036933e-002 

6.26401586397507e-004 

The  Newton-Rhapson  algorithm  is  very  sensitive  to  the  precision  of  the  basis  frequency 
frequency  approximations  since  it  can  easily  settle  on  an  adjacent  peak.  This  was  the  case 
for  the  ISS  data.  The  most  prominent  peaks  could  not  be  identified  by  the  /2  frequencies 
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alone.  The  author  directly  examined  the  spectral  lines  to  improve  the  accuracy  of  the 
approximate  basis  set  with  a  manual  override.  Figure  16  -  Figure  18  show  segments  of 
the  power  spectral  density  plots  used  to  identify  the  frequencies. 
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Power, 


Appendix  A  contains  a  full  accounting  of  the  power  spectral  frequency  combinations 


from  0  to  3  rad/TU. 


Figure  16.  PSD  plot  identifying  the  apsidal  frequency  in  the  z  coordinate. 
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Figure  17.  PSD  plot  identifying  the  nodal  frequency  in  the  x  and  y  coordinates. 


Figure  18.  PSD  plot  identifying  the  anomalistic  frequency  in  the  x,  y,  and  z  coordinates. 


110 


il2  and  n3  were  identified  by  producing  PSD  plots  near  the  expected  peaks  with  up 
to  seven  significant  digits,  was  identified  to  a  higher  precision  from  the  Newton- 
Rhapson  algorithm  which  successfully  identified  the  /  =  (1 , 0 , 1)T  peak  in  previous 
software  outputs.  Deducting  il3  from  this  frequency  provided  a  good  approximation  of 
nx.  Thus,  the  basis  frequencies  that  were  used  to  initialize  the  spectral  analysis  are  given 
in  Table  3. 


Table  3.  Manual  basis  frequencies  for  spectral  analysis. 


Coordinate 

Frequency  (rad/TU) 

9. 234041 79345446e-001 

^2 

-5.967015e-002 

6.273425e-004 

It  is  important  to  note  that  the  manual  basis  set  is  not  the  final  solution;  it  is  merely  an 
approximation  used  to  identify  prominent  spectral  lines.  The  frequencies  from  those  lines 
are  used  to  generate  a  least  squares  estimate  of  the  basis  set.  A  quick  review  of  the 
spectral  lines  employed,  listed  in  Table  4,  shows  that  the  line  peak  powers  are  reasonable 
and  their  frequencies  are  extremely  close  to  the  expected  frequencies  from  the  manual 
set.  The  estimated  basis  set  from  the  least  squares  solution  is  shown  in  Table  5. 

Table  4.  Line  peak  analysis  from  NAFF  and  Newton-Rhapson  results. 


Native 

Coordinate 

Spectral 
Line  (f) 

Residual  Freq.  (rad/TU) 
Actual  -  Expected 

Power  Spectral  Density,  |$|2 

X 

(1,1,1) 

2.04460104402671e-009 

1.82588844891577e-001 

(1,-1, 1) 

-2. 0455929172  7691e-009 

9.99713308702777e-003 

(2,1,1) 

1.844653 11980908e-006 

1.38656419943255e-008 

(2, -1,1) 

1.82975753992842e-006 

7.67087643648475e-010 

y 

(1,1,1) 

2.04458305841371e-009 

1.825853 11258691e-001 

(1,-1,1) 

-2.04558192606896e-009 

9.99362401 708964e-003 

(2,1,1) 

1.84063454988781e-006 

1.38656758742673e-008 

(2, -1,1) 

1.86371049992751e-006 

7.66612586491125e-010 

z 

(1,0,1) 

0 

1.70839929 17232  le-001 

(2,0,1) 

1.84468687014494e-006 

1.30338647593972e-008 

Ill 


Table  5.  Estimated  basis  frequencies. 


Coordinate 

Frequency  (rad/TU) 

9.23406024034362e-001 

^2 

-5.96701449186667e-002 

6.25497810686326e-004 

As  one  final  check,  the  estimated  basis  set  was  used  to  generate  new  expected 
frequencies  for  the  ten  most  prominent  spectral  lines  to  ensure  that  the  residuals  are  small 
after  differencing  the  actual  and  expected  frequencies.  Table  6  depicts  these  values  and, 
again,  they  are  quite  convincing. 

Table  6.  Residuals  from  actual  and  estimated  basis  frequencies. 


Native 

Coordinate 

Spectral  Line  (f) 

Residual  Freq.  (rad/TU) 
Actual  -  Expected 

X 

(1,1,1) 

2.0449984e-009 

(1,-1,1) 

2.0449805e-009 

(2,1,1) 

-2.0451945e-009 

(2, -1,1) 

-2.0451836e-009 

y 

(1,1,1) 

3.9779291e-013 

(1,-1,1) 

-3.5392383e-011 

(2,1,1) 

-4.0539610e-009 

(2, -1,1) 

-1.4930977e-008 

Z 

(1,0,1) 

1.9021981e-008 

(2,0,1) 

-1.6506796e-012 

Given  a  seemingly  good  fit  to  the  spectral  data,  the  estimated  frequencies  were  used 
to  extract  the  Fourier  coefficients.  The  resulting  Fourier  series,  truncated  to  order 
M  =  (6,  17,  6),  was  used  to  generate  native  coordinates  to  compare  to  the  integrated 
orbit.  The  limit  on  the  index  summation  vector  was  decided  after  numerous  trials  with 
lower  and  higher  order  truncations.  The  Fourier  series  from  lower  order  truncations  could 
be  generated  and  evaluated  faster,  but  also  introduced  greater  error.  Craft  shows  the 
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relationship  between  improvements  in  accuracy  and  increases  in  the  number  of  terms 
[52:40].  At  a  certain  point,  the  improvements  become  computationally  burdensome  with 
little  gain  in  accuracy,  if  at  all.  He  also  showed  that  lower  altitude  orbits  require  more 
terms  to  compensate  for  stronger  geopotential  perturbations.  At  an  altitude  of  320  km 
(which  is  close  enough  to  the  ISS  altitude  near  350  km),  he  suggests  that  there  is  little 
relative  accuracy  gain  after  750  frequency  combinations.  A  truncation  of  order  M  = 
(6,  17,  6)  corresponds  to  2,958  frequency  combinations,  far  exceeding  his  accuracy 
threshold  where  errors  are  sub-meter.  Since  the  time  and  effort  had  already  been  invested 
in  generating  this  torus  of  higher  order  truncation,  it  made  little  sense  to  downgrade  to  a 
torus  produced  with  lower  order  terms.  Doing  so  would  have  sacrificed  higher  precision 
for  nearly  unnoticeable  timing  improvements  in  the  MATLAB  code  that  sums  the  Fourier 
series. 

With  the  chosen  truncation,  Figure  19  shows  the  RMS  residuals  for  coordinates  x,  y 
and  z  are  121.9,  121.7,  and  75.7  meters,  respectively.  Since  the  truncation  is  not  a 
significant  cause  of  the  error  and  since  a  periodic  nature  appears  in  the  residuals,  it  must 
be  due  to  incorrect  basis  frequencies. 

Since  the  author  is  confident  that  the  Newton-Rhapson  algorithm  nailed  the  peak  at 
a)i  o,i  =  SI  ■  (1 , 0 , 1 Y ,  fl2  and  ft3  are  the  only  frequencies  in  question  since  is  just  a 
function  of  H3,  specifically  =  o>i0,i  —  ft3.  To  improve  the  frequencies  without 
repeating  the  spectral  analysis,  a  honing  function  was  written  that  adjusted  the 
frequencies  until  the  RMS  residuals  were  minimized.  This  is  a  two-step  process  that  first 
tweaks  fl3  while  holding  £L2  constant  and  allowing  to  vary  with  fi3.  In  the  second 
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Figure  19.  Residuals  from  differencing  the  integrated  orbit  and  reference  torus  native 
coordinates  (uses  Fourier  series  truncation  of  order  M  =  (6,  17,  6)  and  NAFF 

estimated  frequencies). 


step,  Oi  and  fl3  are  held  constant  while  fl2  is  tweaked.  fl2  was  initialized  at  the  NAFF 
estimate  and  fl3  was  initialized  at  the  value  previously  determined  by  detailed  PSD  plots 
near  the  expected  peak.  It  really  doesn’t  matter  where  they  are  initialized,  but  a  smart 
choice  does  speed  up  the  process.  The  honed  basis  set  and  the  residuals  from  the  previous 
estimates  are  given  in  Table  7. 

Table  7.  Honed  basis  frequencies  and  residuals  from  NAFF  estimate. 


Coordinate 

Honed  Basis  Set  (rad/TU) 

Residuals  (rad/TU) 

Honed  -  Estimated 

9.23407979345446  e-001 

1.95531 1083046 14e-006 

^2 

-5.96701428186668  e-002 

2. 100000 19457146e-009 

fi3 

6.235425  e-004 

-1.9553 1068455002e-006 
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The  newly  honed  basis  set  was  used  to  test  the  deterministic  torus  accuracy  again. 
Figure  20  shows  the  RMS  residuals  for  coordinates  x,  y  and  z  are  18.7,  17.6,  and  19.8 
meters,  respectively.  Although  the  result  is  a  significant  improvement  over  the  previously 
estimated  frequencies,  it  must  be  noted  that  the  coefficients  still  contain  errors  from 
picking  off  peaks  using  the  NAFF  estimated  frequencies.  A  better  fit  between  the  torus 
and  the  integrated  orbit  cannot  be  obtained  until  the  NAFF  method  is  improved  to 
identify  a  more  accurate  basis  set. 


Figure  20.  Residuals  from  differencing  the  integrated  orbit  and  reference  torus  native 
coordinates  (uses  Fourier  series  truncation  of  order  M  =  (6,  17,  6)  and  honed 

frequencies). 

As  a  final  comment,  one  interesting  observation  from  a  review  of  the  spectra  in  the 
frequency  domain  (see 


115 


Appendix  A)  is  that  the  peaks  used  for  the  least  squares  solution  given  in  Table  4  are 
not  the  most  prominent  peaks.  The  most  prominent  peaks  are  centered  about  o)n  0  n  =  ft  ■ 
(n,0  ,n)T  in  the  z  coordinate.  For  example,  the  peak  thought  to  be  the  maximum  at 
(o201  =  ft  ■  (2 , 0 , 1)T  is  actually  overshadowed  by  a>2 0,2  =  ft  ■  (2 , 0  , 2)r.  Similarly, 
the  x  and  y  coordinates  show  the  same  behavior.  The  peak  thought  to  be  the  maximum  at 
<^2,1,1  =  ft  ■  (2 , 1 , 1  )r  is  actually  overshadowed  by  (o212  =  ft  ■  (2 , 1 , 2)r.  This  may  be 
the  reason  the  least  squares  solution  gives  the  wrong  frequencies. 

4.2  Filter  Verification 

A  verification  of  the  NLS  and  Bayes  filters  was  completed  using  observation  data  that 
was  manufactured  directly  from  the  reference  torus  using  a  constant  SY0  offset.  This  was 
necessary  to  show  that  the  filters  converge  on  the  correct  <5T0  and  can  be  trusted  with  real 

data  when  £T0  is  unknown.  The  data,  shown  in  Figure  21,  was  generated  from  the 

following  offset  to  the  reference  torus: 

0.003051846308510 
3.364795882078831  x  10"6 
-0.003041378665051 
2.530687412213637  X  10-6 
-2.032460083183908  x  10"5 
2.521736621432707  X  10-6 

The  known  offsets  were  recovered  by  the  filters  to  within  machine  double  precision. 
As  confirmation  that  the  filter  indeed  found  the  correct  offsets  to  the  reference  torus, 
Figure  22  shows  the  native  coordinate  residuals  from  the  manufactured  observation  data 

and  the  updated  reference  torus  (updated  with  the  discovered  <SK0). 
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Figure  21.  Residuals  from  manufactured  observation  data  generated  by  offsets  to  the 

reference  torus. 
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Figure  22.  Coordinate  residuals  show  reference  torus  was  corrected  to  match  the 

manufactured  data. 
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Phase  Angle  (rad) 


Of  further  interest  is  the  effect  the  simulated  offsets  had  on  the  reference  torus.  The 


filter  results  show  a  linear  drift  in  the  initial  phase  angles  from  the  SQ  updates  at  each 
new  epoch  as  indicated  in  Figure  23  -  Figure  25.  This  is  expected  from  (172). 


Time  (days) 

Figure  23.  Linear  drift  in  QOi  due  to  constant  dPi. 
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Figure  24.  Linear  drift  in  QO2  due  to  constant  dP2. 
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Figure  25.  Linear  drift  in  QO3  due  to  constant  dP3. 
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4.3  Orbit  Fitting 

The  KAM  torus  was  used  to  generate  deterministic  and  stochastic  predictions  to 
compare  to  the  ISS  data.  The  predictions  from  the  purely  deterministic  model  gives  a  real 
sense  of  how  quickly  the  ISS  navigates  off  of  the  reference  torus  without  perturbation 
updates.  The  stochastic  predictions  were  formed  with  updates  to  the  reference  torus  using 
a  pseudo-stochastic  parameterization  of  the  initial  phase  angles  and  momenta  offsets. 

4.3.1  Deterministic  Predictions 

In  its  unmodified  state  at  epoch,  the  reference  torus  can  be  used  to  generate 
deterministic  predictions  of  the  ISS  native  coordinates  when  time  is  the  only  parameter 
that  is  varied.  This  allows  for  a  direct  comparison  to  the  actual  ISS  data  to  discern  how 
well  the  gravitationally  perturbed  torus  represents  reality.  Doing  so  reveals  a  quadratic 
growth  in  the  residuals  exceeding  500  km  in  just  8  days.  Figure  26  shows  the  residuals 
from  the  raw  data,  whereas  Figure  27  shows  the  residuals  after  fallacious  data  is  rejected 
by  the  pre-filter.  Notice  that  the  pre-filter  removed  a  lot  of  trashy  data  near  the  11th  day. 
This  is  an  important  clue  for  later  analysis  of  the  stochastic  predictions. 
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Figure  26.  Residuals  from  direct  comparison  of  deterministic  reference  torus  to 

observation  data. 


Figure  27.  Pre-filtered  residuals  from  direct  comparison  of  deterministic  reference  torus 

to  observation  data. 
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4.3.2  Stochastic  Predictions 


The  time  varying  polynomials  characterizing  the  initial  phase  angles  and  momenta 
offsets  were  used  to  generate  stochastic  predictions  from  the  reference  torus.  Since  their 
behavior  may  be  sensitive  to  the  batch  size  during  the  estimation  process,  batches  of  100, 
300  and  1000  observations  were  used  at  roughly  1  minute  intervals.  Only  the  first  half  of 
the  ISS  data  were  used  in  the  fitting  process,  whereas  the  remaining  half  were  used  as 
fresh,  new  observations  to  simulate  a  “real  time”  evaluation  of  the  predictions. 

4.3.2.1  Estimations  and  Predictions:  Batches  of  1000  Observations 

As  an  assurance  that  the  data  were  processed  correctly  by  the  pre-filter,  NLS  and 
Bayes  filter,  results  from  the  first  two  batches  will  be  presented  in  addition  to  the  full 
string  of  21  batches  required  to  process  the  nearly  2-weeks  of  ISS  data  at  1,000 
observations  per  batch. 

The  first  batch  of  data  was  processed  by  the  NLS  since  a  priori  data  was  not 
available.  Prior  to  the  first  iteration  of  NLS,  the  residuals  from  the  raw  data  are  plotted  in 
Figure  28  for  monitoring  the  accuracy  of  the  pre-filter.  Figure  29  shows  the  residuals 
after  pre-filtration  and  Figure  30  shows  the  residuals  after  10  iterations  of  NLS.  The  final 
RMS  residuals  for  coordinates  x,  y  and  z  are  357.17,  334.23,  and  338.26  meters, 
respectively. 
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gure  28.  Residuals  from  first  set  of  1000  observations  prior  to  pre-filter  and  NLS. 


Figure  29.  Pre-filtered  residuals  from  first  set  of  1000  observations  prior  to  NLS. 
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Figure  30.  NLS  filtered  residuals  from  first  set  of  1000  observations. 


Subsequent  batches  (after  the  NLS  batch)  are  processed  by  the  Bayes  filter.  The 
residuals  from  the  second  batch  of  raw  data  are  plotted  in  Figure  31  for  monitoring  the 
accuracy  of  the  pre-filter  prior  to  the  first  iteration  of  the  Bayes  filter.  Figure  32  shows 
the  residuals  after  pre-filtration  and  Figure  33  shows  the  residuals  after  10  iterations  of 
Bayes.  The  final  RMS  residuals  for  coordinates  x,  y  and  z  are  471.53,  326.24,  and  263.47 
meters,  respectively. 
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Figure  31.  Residuals  from  second  set  of  1000  observations  prior  to  pre -filter  and  Bayes. 
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Figure  32.  Pre-filtered  residuals  from  second  set  of  1000  observations  prior  to  Bayes. 
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Figure  33.  Bayes  filtered  residuals  from  second  set  of  1000  observations. 

A  concatenation  of  the  final  residuals  after  10  iterations  from  all  21  batches  is  shown 
in  Figure  34.  The  RMS  residuals  for  coordinates  x,  y  and  z  are  495.5,  440.9,  and  370.1 
meters,  respectively.  The  same  is  done  for  the  reference  torus  corrections  in  Figure  35  - 
Figure  37  to  show  the  change  in  the  initial  phase  angles  and  momenta  that  best  fit  the 
data.  The  plots  reveal  a  linear  trend  in  the  momenta  and  a  periodicity  in  the  phase  angles. 
The  linearity  in  the  momenta  is  a  pleasant  finding  since  a  line  can  easily  be  fit  to  the  data 
-  precisely  what  is  needed  to  utilize  the  reference  torus  for  stochastic  predictions.  The 
periodicity  in  the  phase  angle  corrections  is  not  cause  for  alarm  since  they  oscillate  about 
zero  with  small  amplitudes  which  will  have  a  marginal  effect  when  pooled  with  the  initial 
phase  angles,  but  a  peculiar  pattern  arises  in  which  dQi  and  dQ3  exactly  oppose  one 
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another.  As  dQi  increases,  dQ3  decreases,  and  vice  versa.  This  is  the  classic  mean 
anomaly  vs.  argument  of  perigee  problem  for  low  eccentricity  orbits. 

Recall  that  the  time  derivative  of  Qi  is  analogous  to  mean  motion  and  the  time 
derivative  of  Q3  is  analogous  to  the  apsidal  regression  rate,  so  both  coordinates  have  the 
physical  effect  of  displacements  on  the  orbital  plane.  Since  the  ISS  orbit  has  a  very  low 
eccentricity  (-0.0007  -  0.002)  the  apsidal  regression  rate  is  very  small  (~7.7e-7  s'1)  and  is 
easily  overlooked  by  the  mean  motion  (-0.00114  s'1)  in  the  same  plane.  This  makes  it 
susceptible  to  modeling  errors  as  the  two  coordinates  fight  to  converge  upon  a  solution. 
More  will  be  said  on  this  later. 


Figure  34.  Bayes  filtered  residuals  from  batches  of  1000  observations. 
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Figure  35.  dQi  and  dPi  torus  corrections  from  Bayes  filtered  batches  of  1000 

observations. 
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Figure  36.  dQ2  and  dP2  torus  corrections  from  Bayes  filtered  batches  of  1000 

observations. 
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Figure  37.  dQ3  and  dP3  torus  corrections  from  Bayes  filtered  batches  of  1000 

observations. 


After  accumulating  the  updates  to  Q0  from  (172),  a  polynomial  of  degree  two  was  fit 

to  the  data.  Similarly,  a  polynomial  of  degree  one  was  fit  to  SP  as  shown  in  Figure  38  - 
Figure  40. 
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Figure  38.  Polynomial  approximations  of  QOi  and  dPi  from  Bayes  filtered  batches  of 

1000  observations. 
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Figure  39.  Polynomial  approximations  of  QO2  and  dP2  from  Bayes  filtered  batches  of 

1000  observations. 
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Figure  40.  Polynomial  approximations  of  QO3  and  dP3  from  Bayes  filtered  batches  of 

1000  observations. 


The  polynomial  functions  representing  the  pseudo-stochastic  parameters  were  used 
with  the  reference  torus  to  generate  stochastic  predictions.  All  42,000  observations, 
including  the  first  21,000  that  were  used  during  the  estimation  process,  were  compared  to 
the  stochastic  prediction  from  the  torus.  The  last  21,000  observations  are  used  to  simulate 
the  torus’  ability  to  predict  in  “real  time”  since  these  data  were  not  processed  previously 
by  the  Bayes  filter.  Figure  41  shows  the  residuals  for  the  full  42,000  observations.  Prior 
to  the  degradation  in  the  prediction  at  18.5  days,  the  RMS  residuals  for  coordinates  x,  y 
and  z  are  2.586,  2.602,  and  2.455  kilometers,  respectively. 
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Figure  41.  Residuals  generated  using  pseudo-stochastic  parameters  estimated  from 

batches  of  1000  observations. 


After  observing  the  results,  the  obvious  question  is:  why  did  the  predictions  hold  so 
consistently  for  18.5  days  (almost  4  days  into  the  real  time  segment)  and  then  suddenly 
breakdown  midway  through  April  12?  A  review  of  the  dynamic  events  during  the 
prediction  times,  listed  in  Table  8,  show  no  activity  on  that  day.  However,  a  correlation 
of  these  events  with  the  results  gives  an  important  clue.  Consider  again  the  residuals 
depicted  with  these  overlapping  events  in  Figure  42. 

Table  8.  ISS  Dynamic  Events. 


Time  (GMT) 

Event 

4  Apr  03:10-06:50 

Soyuz  (TMA-18)  docking 

7  Apr  04:13-07:44 

Discovery  (STS- 131)  docking 

17  Apr  10:20-15:06 

Discovery  (STS-131)  undocking 

17  Apr  19:55  -  18  Apr  02:15 

Progress  35  Prop  Purge 

22  Apr  13:15-18:04 

Progress  35  undocking 
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Figure  42.  Residuals  from  stochastic  predictions  correlated  with  dynamic  events. 


velocity  vector  (XVV  attitude) 
and  the  z-axis  pointing  nadir  [53], 
This  means  the  velocity  vector  is 
nearly  normal  to  the  underside  of 
the  shuttle.  Figure  43  gives  a  sense 
of  how  much  additional  drag  is 


Notice  that  the  slope  of  the  residuals  changes  at  the  time  of  the  Shuttle’s  undocking. 
This  is  due  to  a  large  decrease  in  drag  that  causes  a  divergence  from  the  pseudo¬ 
stochastic  torus  parameters.  In  its 
current  configuration,  the  ISS  is 
flown  with  its  x-axis  in  the 


Figure  43.  Artist  rendering  of  the  Space  Shuttle 
docked  with  the  ISS.  Credit:  NASA 
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incurred  when  the  shuttle  is  docked. 


It  stands  to  reason  that  any  significant  change  in  drag  or  increase  in  thrust  will  cause  a 
divergence  from  the  predictions.  The  arrival  of  the  Soyuz  and  Shuttle  did  not  induce  the 
same  divergence  because  the  pseudo-stochastic  parameters  were  estimated  during  those 
docking  events  and  therefore  the  curves  were  bent  to  compensate  for  the  change  in  drag. 
But  this  does  not  excuse  the  Russian  Progress  cargo  ship’s  departure  from  the  ISS  on 

April  22.  Figure  44  shows  a 
graphic  of  the  ISS  with  the 
solar  arrays  and  Active 
Thermal  Control  System 
(ATCS)  radiators  removed 
for  viewing  of  the  Russian 
segment,  particularly  a 
Russian  Progress  resupply 
vehicle  docked  at  the  Poisk 
Mini  Research  Module  2  (MRM2).  This  zenith  port  was  the  same  one  used  for  Progress 
35  so  it  gives  a  sense  of  the  drag  contribution.  Because  its  projected  surface  area  is  so 
much  smaller  than  the  Shuttle,  its  effect  is  more  difficult  to  observe  in  the  residual  growth 
that  is  already  increasing  rapidly  due  to  the  Shuttle’s  departure. 

The  Shuttle’s  departure  is  proof  that  a  significant  change  in  drag  will  cause 
degradation  in  the  torus  prediction,  but  if  there  were  no  changes  to  the  total  surface  area 
of  the  ISS  on  April  12,  the  only  other  possibility  for  a  change  in  drag  is  a  change  in 


Figure  44.  Scaled  comparison  of  Progress  resupply  vehicle 
and  Space  Shuttle  docked  to  the  ISS.  Credit:  NASA 
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particle  density.  Space  weather  reports  from  the  National  Oceanic  and  Atmospheric 


Administration  (NOAA)  were  used  to  assess  geomagnetic  activity  that  may  have 
contributed  to  changes  in  drag  that  day  [54], 

The  NOAA  uses  the  K-index  from  the  Boulder  Magnetometer  to  approximate  the 
planetary  Kp-index  for  real-time  alerts  and  warnings.  According  to  the  NOAA  Space 
Weather  Scale  for  Geomagnetic  Storms  [55],  minor  storms  (G1  class)  are  characterized 
by  Kp  values  less  than  or  equal  to  5.  These  typically  have  negligible  impacts  on 
spacecraft  operations.  Moderate  storms  (G2  class)  are  characterized  by  Kp  values  equal 
to  6.  Under  these  conditions  the  NOAA  states,  “corrective  actions  to  orientation  may  be 
required  by  ground  control;  possible  changes  in  drag  affect  orbit  predictions”  [55].  Strong 
storms  (G3  class),  severe  storms  (G4  class)  and  extreme  storms  (G5  class)  are 
characterized  by  Kp  values  equal  to  7,  8  (including  a  9-)  and  9,  respectively.  The  NOAA 
warns  that  the  effects  of  these  higher  class  storms  usually  include  surface  charging, 
orientation,  tracking  and  prediction  problems  commensurate  with  the  severity  of  the 
storm.  G3  class  storms  are  powerful  enough  that  the  NOAA  specifically  warns  LEO 
satellites  can  experience  increased  drag.  For  the  time  period  of  ISS  observations,  all 
NOAA  alerts  for  G2  through  G5  class  storms  were  logged  in  Table  9  [56-60]. 

Table  9.  Moderate  to  Extreme  Space  Weather  Events. 


1  Time  (GMT) 

Event 

05  Apr  0920 

ALERT:  Geomagnetic  K  =  6 

05  Apr  0955 

ALERT:  Geomagnetic  K  =  7 

06  Apr  0422 

ALERT:  Geomagnetic  K  =  6 

12  Apr  0225 

ALERT:  Geomagnetic  K  =  6 

12  Apr  0240 

ALERT:  Geomagnetic  K  =  7 
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From  25  March  through  4  April 
2010,  the  NOAA  reported  mostly 
quiet  to  minor  storm  levels  in  the 
geomagnetic  field  with  only  brief 
active  periods  where  the  Kp  index 
peaked  at  4.  A  halo  coronal  mass 
ejection  (CME)  observed  on  3 
April  caused  an  increase  in 
geomagnetic  activity  starting  on  5 

April  that  peaked  with  G3  class 
storms.  The  effects  of  the  CME 
waned  on  7  April.  Any  change  in 
particle  density  incurred  from  this  CME  would  have  been  absorbed  by  the  pseudo¬ 
stochastic  parameters  during  the  estimation  process  and  won’t  appear  in  the  predictions. 
Conditions  were  then  mostly  quiet  with  some  isolated  minor  storm  ensuing  through  11 
April.  A  full  halo  CME  aimed  almost  directly  at  earth  was  observed  on  8  April  at  0325 
UTC  causing  major  to  sever  storms  on  12  April.  It  is  shown  in  Figure  45.  By  13  April  the 
active  storms  dwindled  to  quiet  levels  at  which  conditions  predominantly  remained  with 
only  minor  storms  through  24  April. 

The  presence  of  major  geomagnetic  storms  on  12  April  certainly  leaves  a  possible 
explanation  for  the  degradation  in  the  predictions,  but  to  be  sure,  further  investigations 
will  be  required.  The  first  step  is  to  run  the  last  set  of  21,000  observations  through  the 


Figure  45.  Solar  flare  in  upper  left  quadrant  on 
8  April  2010  produced  an  earth-directed  CME. 
Credit:  NASA  Solar  Dynamics  Observatory 
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Bayes  filter  to  estimate  what  the  pseudo-stochastic  parameters  should  have  been  and 
compare  those  to  the  curves  generated  from  the  first  set  of  21,000  observations. 

Figure  46  -  Figure  48  show  that  the  polynomials  fit  the  estimates  fairly  well  until 
April  12  when  the  predictions  degenerated.  A  clear  jump  is  observed  in  the  momenta 
offsets  at  the  times  of  interest,  further  supporting  the  hypothesis  that  drag  changes  shifted 
the  ISS  onto  an  adjacent  torus  trajectory  that  was  not  modeled. 
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Figure  46.  Estimates  of  dPi  from  Bayes  filtered  batches  of  1000  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  47.  Estimates  of  dP2  from  Bayes  filtered  batches  of  1000  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  48.  Estimates  of  CIP3  from  Bayes  filtered  batches  of  1000  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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The  indisputable  proof  is  found  by  calling  upon  the  unforsaken  osculating  orbital 
elements.  The  semimajor  axis  is  known  to  transfigure  with  variations  in  atmospheric 
drag.  Figure  49  shows  that  this  is  precisely  what  happened.  Because  the  change  is  so 
subtle,  the  event  markups  have  been  removed  in  Figure  50  to  more  clearly  reveal  the 
changes. 
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Figure  49.  The  pre-filtered,  osculating  semimajor  axis  reveals  subtle  mutations  near  the 

events  of  interest. 
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Figure  50.  The  pre-filtered,  osculating  semimajor  axis  reveals  subtle  mutations  near  the 
events  of  interest.  The  event  markups  have  been  removed  for  unaided  viewing. 

The  next  set  of  results  will  reveal,  prior  to  the  time  at  which  the  predictions  break 
down  due  to  the  G3-class  storm,  whether  or  not  the  predictions  can  be  improved  using 
pseudo-stochastic  parameters  estimated  with  different  sized  batches. 

4.3.2.2  Estimations  and  Predictions:  Batches  of  300  Observations 

The  ensuing  results  are  produced  from  batches  of  300  observations  which  translate  to 
a  little  more  than  3  full  orbital  revolutions.  As  an  assurance  that  the  data  was  processed 
correctly  by  the  pre-filter,  NLS  and  Bayes  filter,  results  from  the  first  two  batches  will  be 
presented  in  addition  to  the  full  string  of  70  batches  required  to  process  the  nearly  2- 
weeks  of  ISS  data  at  300  observations  per  batch. 
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The  first  batch  of  data  was  processed  by  the  NLS  since  a  priori  data  was  not 
available.  Prior  to  the  first  iteration  of  NLS,  the  residuals  from  the  raw  data  are  plotted  in 
Figure  51  for  monitoring  the  accuracy  of  the  pre-filter.  Figure  52  shows  the  residuals 
after  pre-filtration  and  Figure  53  shows  the  residuals  after  10  iterations  of  NLS.  The  final 
RMS  residuals  for  coordinates  x,  y  and  z  are  225.49,  202.20,  and  195.93  meters, 
respectively. 
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Figure  51.  Residuals  from  first  set  of  300  observations  prior  to  pre-filter  and  NLS. 
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Figure  52.  Pre-filtered  residuals  from  first  set  of  300  observations  prior  to  NLS. 


Figure  53.  NLS  filtered  residuals  from  first  set  of  300  observations. 
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Subsequent  batches  (after  the  NLS  batch)  are  processed  by  the  Bayes  filter.  The 
residuals  from  the  second  batch  of  raw  data  are  plotted  in  Figure  54  for  monitoring  the 
accuracy  of  the  pre-filter  prior  to  the  first  iteration  of  the  Bayes  filter.  Figure  55  shows 
the  residuals  after  pre-filtration  and  Figure  56  shows  the  residuals  after  10  iterations  of 
NLS.  The  final  RMS  residuals  for  coordinates  x,  y  and  z  are  218.91,  118.23,  and  121.68 
meters,  respectively. 


Figure  54.  Residuals  from  second  set  of  300  observations  prior  to  pre-filter  and  Bayes. 
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gure  55.  Pre-filtered  residuals  from  second  set  of  300  observations  prior  to  Bayes. 
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Figure  56.  Bayes  filtered  residuals  from  second  set  of  300  observations. 


144 


A  concatenation  of  the  final  residuals  after  10  iterations  from  all  70  batches  is  shown 


in  Figure  57.  The  RMS  residuals  for  coordinates  x,  y  and  z  are  255.14,  243.40,  and 
224.80  meters,  respectively.  The  same  is  done  for  the  reference  torus  corrections  in 
Figure  58  -  Figure  60  to  show  the  change  in  the  initial  phase  angles  and  momenta  that 
best  fit  the  data.  The  plots  reveal  a  linear  trend  in  the  momenta  and  a  periodicity  in  the 
phase  angles  as  observed  and  discussed  previously  in  §4.3.2. 1. 


Figure  57.  Bayes  filtered  residuals  from  batches  of  300  observations. 
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Figure  58.  dQi  and  dPi  torus  corrections  from  Bayes  filtered  batches  of  300 

observations. 
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Figure  59.  dQ2  and  dP2  torus  corrections  from  Bayes  filtered  batches  of  300 

observations. 
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Figure  60.  dQ3  and  dP3  torus  corrections  from  Bayes  filtered  batches  of  300 

observations. 


After  accumulating  the  updates  to  Q0  from  (172),  a  polynomial  of  degree  two  was  fit 

to  the  data.  Similarly,  a  polynomial  of  degree  one  was  fit  to  SP  as  shown  in  Figure  61  - 
Figure  63. 
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Figure  61.  Polynomial  approximations  of  QOi  and  dPi  from  Bayes  filtered  batches  of  300 

observations. 
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Figure  62.  Polynomial  approximations  of  QO2  and  dP2  from  Bayes  filtered  batches  of  300 

observations. 
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Figure  63.  Polynomial  approximations  of  QO3  and  dP3  from  Bayes  filtered  batches  of  300 

observations. 


As  before,  the  pseudo-stochastic  parameters  were  used  to  generate  stochastic 
predictions  from  the  reference  torus.  All  42,000  observations,  including  the  first  21,000 
that  were  used  during  the  estimation  process,  were  compared  to  the  torus  prediction.  The 
last  21,000  observations  are  used  to  simulate  the  torus’  ability  to  predict  in  real  time  since 
the  data  were  not  processed  previously  by  the  Bayes  filter.  Figure  64  shows  the  residuals 
for  the  full  42,000  observations  with  the  exact  same  deterioration  characteristics  as  seen 
previously.  Prior  to  the  degradation  in  the  prediction  at  18.5  days,  the  RMS  residuals  for 
coordinates  x,  y  and  z  are  2.606,  2.617,  and  2.456  kilometers,  respectively.  These  are 
only  slightly  worse  than  those  found  with  batches  of  1,000  observations. 
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Figure  64.  Residuals  generated  using  pseudo-stochastic  parameters  estimated  from 

batches  of  300  observations. 


For  thoroughness,  the  last  set  of  21,000  observations  were  also  run  through  the  Bayes 
filter  to  estimate  what  the  pseudo-stochastic  parameters  should  have  been  and  compare 
those  to  the  curves  generated  from  the  first  set  of  21,000  observations.  In  doing  so,  the 
same  trends  appear  as  before  with  1,000  observations  per  batch.  Figure  65  -  Figure  67 
show  that  the  polynomials  fit  the  estimates  fairly  well  until  April  12  when  the  predictions 
degenerated.  A  clear  jump  is  observed  in  the  momenta  offsets  at  the  times  of  interest, 
indicating  that  the  batch  size  likely  has  nothing  to  do  with  the  anomalies.  This  further 
supports  the  hypothesis  that  drag  changes  shifted  the  ISS  onto  an  adjacent  torus  trajectory 
that  was  not  modeled. 
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Figure  65.  Estimates  of  dPi  from  Bayes  filtered  batches  of  300  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  66.  Estimates  of  dP2  from  Bayes  filtered  batches  of  300  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  67.  Estimates  of  dP3  from  Bayes  filtered  batches  of  300  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 


4.3.2.3  Estimations  and  Predictions:  Batches  of  100  Observations 

The  ensuing  results  are  produced  from  batches  of  100  observations  which  translate  to 
a  little  more  than  a  full  orbit.  As  an  assurance  that  the  data  was  processed  correctly  by  the 
pre-filter,  NLS  and  Bayes  filter,  results  from  the  first  two  batches  will  be  presented  in 
addition  to  the  full  string  of  210  batches  required  to  process  the  nearly  2-weeks  of  ISS 
data  at  100  observations  per  batch. 

The  first  batch  of  data  was  processed  by  the  NLS  since  a  priori  data  was  not 
available.  Prior  to  the  first  iteration  of  NLS,  the  residuals  from  the  raw  data  are  plotted  in 
Figure  68  for  monitoring  the  accuracy  of  the  pre-filter.  Figure  69  shows  the  residuals 
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after  pre-filtration  and  Figure  70  shows  the  residuals  after  10  iterations  of  NLS.  The  final 


RMS  residuals  for  coordinates  x,  y  and  z  are  104.95,  51.44,  and  99.25  meters, 
respectively. 
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Figure  68.  Residuals  from  first  set  of  100  observations  prior  to  pre-filter  and  NLS. 
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igure  69.  Pre-filtered  residuals  from  first  set  of  100  observations  prior  to  NLS. 


Figure  70.  NLS  filtered  residuals  from  first  set  of  100  observations. 
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Subsequent  batches  (after  the  NLS  batch)  are  processed  by  the  Bayes  filter.  The 
residuals  from  the  second  batch  of  raw  data  are  plotted  in  Figure  71  for  monitoring  the 
accuracy  of  the  pre-filter  prior  to  the  first  iteration  of  the  Bayes  filter.  Figure  72  shows 
the  residuals  after  pre-filtration  and  Figure  73  shows  the  residuals  after  10  iterations  of 
NLS.  The  final  RMS  residuals  for  coordinates  x,  y  and  z  are  80.63,  68.75,  and  97.30 
meters,  respectively. 


Figure  71.  Residuals  from  second  set  of  100  observations  prior  to  pre-filter  and  Bayes. 
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ure  72.  Pre-filtered  residuals  from  second  set  of  100  observations  prior  to  Bayes. 


Figure  73.  Bayes  filtered  residuals  from  second  set  of  100  observations. 
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A  concatenation  of  the  final  residuals  after  10  iterations  from  all  210  batches  is  shown 


in  Figure  74.  The  RMS  residuals  for  coordinates  x,  y  and  z  are  180.42,  181.70,  and 
169.74  meters,  respectively.  The  same  is  done  for  the  reference  torus  corrections  in 
Figure  75  -  Figure  77  to  show  the  change  in  the  initial  phase  angles  and  momenta  that 
best  fit  the  data.  The  plots  reveal  a  linear  trend  in  the  momenta  and  a  periodicity  in  the 
phase  angles  as  observed  and  discussed  previously  in  §4.3.2. 1. 
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Figure  74.  Bayes  filtered  residuals  from  batches  of  100  observations. 
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Figure  75.  dQi  and  dPi  torus  corrections  from  Bayes  filtered  batches  of  100 

observations. 
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Figure  76.  dCh  and  dP2  torus  corrections  from  Bayes  filtered  batches  of  100 

observations. 
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Figure  77.  dQ3  and  dP3  torus  corrections  from  Bayes  filtered  batches  of  100 

observations. 


After  accumulating  the  updates  to  Q0  from  (172),  a  polynomial  of  degree  two  was  fit 

to  the  data.  Similarly,  a  polynomial  of  degree  one  was  fit  to  SP  as  shown  in  Figure  78  - 
Figure  80. 
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Figure  78.  Polynomial  approximations  of  QOi  and  dPi  from  Bayes  filtered  batches  of  100 

observations. 
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Figure  79.  Polynomial  approximations  of  QO2  and  dP2  from  Bayes  filtered  batches  of  100 

observations. 
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Figure  80.  Polynomial  approximations  of  QO3  and  dP3  from  Bayes  filtered  batches  of  100 

observations. 


As  before,  the  pseudo-stochastic  parameters  were  used  to  generate  stochastic 
predictions  from  the  reference  torus.  All  42,000  observations,  including  the  first  21,000 
that  were  used  during  the  estimation  process,  were  compared  to  the  torus  prediction.  The 
last  21,000  observations  are  used  to  simulate  the  torus’  ability  to  predict  in  real  time  since 
the  data  were  not  processed  previously  by  the  Bayes  filter.  Figure  81  shows  the  residuals 
for  the  full  42,000  observations  with  the  exact  same  deterioration  characteristics  as  seen 
previously.  Prior  to  the  degradation  in  the  prediction  at  18.5  days,  the  RMS  residuals  for 
coordinates  x,  y  and  z  are  2.789,  2.800,  and  2.643  kilometers,  respectively.  These  are 
slightly  worse  than  those  found  with  batches  of  300  and  1,000  observations. 
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Figure  81.  Residuals  generated  using  pseudo-stochastic  parameters  estimated  from 

batches  of  100  observations. 


For  thoroughness,  the  last  set  of  21,000  observations  were  also  run  through  the  Bayes 
filter  to  estimate  what  the  pseudo-stochastic  parameters  should  have  been  and  compare 
those  to  the  curves  generated  from  the  first  set  of  21,000  observations.  In  doing  so,  the 
same  trends  appear  as  before  with  300  and  1,000  observations  per  batch.  Figure  82  - 
Figure  84  show  that  the  polynomials  fit  the  estimates  fairly  well  until  April  12  when  the 
predictions  degenerated.  A  clear  jump  is  observed  in  the  momenta  offsets  at  the  times  of 
interest,  further  supporting  the  hypothesis  that  drag  changes  shifted  the  ISS  onto  an 
adjacent  torus  trajectory  that  was  not  modeled. 
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Figure  82.  Estimates  of  dPi  from  Bayes  filtered  batches  of  100  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  83.  Estimates  of  dP2  from  Bayes  filtered  batches  of  100  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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Figure  84.  Estimates  of  CIP3  from  Bayes  filtered  batches  of  100  observations  show  a 
slight  offset  from  the  polynomial  approximations  at  times  of  interest. 
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V.  Conclusions  and  Recommendations 


The  underlying  motivation  for  this  research  was  to  improve  space  situational 
awareness  by  extending  KAM  theorem’s  efficacy  to  stochastic  orbit  prediction.  The 
results  are  conclusive.  Most  stochastic  effects  in  LEO  can  be  modeled  to  fit  and  predict  a 
satellite’s  motion  near  a  reference  torus.  The  geometric  construal  of  the  satellite  path 
suggests  a  flexing  annulus  of  concentric,  continuous  tori  about  which  the  satellite  may 
circumnavigate  systematically. 

5.1  Torus  Construction 

The  Fourier  decomposition  of  the  integrated  trajectory  using  Laskar’s  NAFF  worked 
reasonably  well  to  map  (x,v)  ->  (Q ,P);  however,  significant  issues  must  be  overcome  to 
perfectly  identify  the  fundamental  frequencies.  Even  a  faintly  flawed  basis  set  is  enough 
to  distort  the  Laskar  process  when  picking  off  spectral  amplitudes.  This  is  believed  to 
have  been  the  cause  of  the  inherent  errors  in  the  ISS  reference  that  grew  periodically  to 
roughly  80  meters  within  the  first  six  months. 

As  mentioned  in  §4.1,  the  most  likely  reason  for  the  wrong  frequencies  is  that  the  ten 
most  prominent  spectral  lines  used  for  the  least  squares  solution  were  associated  with  the 

wrong  j-index  label.  The  most  prominent  peaks  are  not  centered  about  con  01  =  ft  ■ 

(n ,  0 , 1)T  in  the  z  coordinate,  but  instead  are  centered  about  <un  0  n  =  ft  ■  (n ,  0  ,  n)T. 
Nevertheless,  these  inconsistencies  in  the  torus-construction  machinery  still  provided  a 
reference  torus  that  was  good  enough  to  proceed  with  the  stochastic  estimation  process  - 
a  testament  of  KAM  theory’s  tractability  even  without  an  exact  basis  set. 
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If  KAM  theory  is  to  find  a  place  in  the  mainstream  of  orbital  sciences,  the  procedures 
for  generating  tori  must  be  improved  so  that  they  may  be  autonomously  generated  from 
space  surveillance  data  with  limited  human  involvement.  The  current  state  of  the  art 
requires  extensive  human  massaging  and  tweaking  that  would  make  it  impractical  as  an 
application  for  hundreds  of  thousands  of  earth  orbiting  objects.  Nevertheless,  it  is  the 
belief  of  this  author  along  with  all  those  who  toil  in  this  field,  that  it  is  possible,  perhaps 
within  this  decade. 

5.2  Stochastic  Predictions 

In  this  present  work  as  in  that  by  Little,  low  earth  orbits  have  exhibited  perturbations 
that  cannot  be  modeled  by  a  purely  deterministic,  static  torus.  If  a  non-chaotic  earth  orbit 
is  occupied,  the  phase  space  will  be  filled  with  a  dense  continuum  of  flexing,  persisting 
tori.  The  presence  of  large  perturbations  tends  to  distort  local  tori,  forcing  the  satellite  to 
migrate  from  one  torus  to  another.  Even  still,  the  results  from  §4.3  showed  that  it  is 
possible  to  forecast  motion  in  the  vicinity  of  a  reference  torus  when  the  right  conditions 
are  met. 

Modeling  the  non-periodic  motion  about  the  reference  torus  is  largely  contingent 
upon  the  pseudo-stochastic  parameters.  The  low  eccentricity  of  the  ISS  made  it  difficult 
to  estimate  the  torus  coordinate  offsets  since  the  Keplarian  and  apsidal  coordinates 
opposed  one  another  in  a  sinusoidal  manner.  Curve-fits  of  the  parameters  were  still 
possible,  but  they  suffered  from  immediate  limitations  in  accuracy.  The  batch  size  had  a 
slight  impact  on  the  accuracy;  larger  batches  provide  better  results,  but  it  is  not  clear  to 
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what  extent  size  matters.  A  more  thorough  review  would  be  required  with  varied 
observation  intervals. 

Another  cause  for  errors  in  the  non-periodic  motion  is  a  rapid  shift  in  the  drag 
environment  or  surface  area.  These  changes  require  larger  offsets  in  the  pseudo-stochastic 
parameters  than  are  modeled  during  the  estimation  process. 

In  the  end,  the  difficulties  encountered  are  minor  compared  to  the  incredible  feat  of 
modeling  the  largest  earth  orbiting  satellite  in  excess  of  18  days  with  RMS  residuals 
bounded  near  2  km  -  quite  the  contrast  from  the  results  that  showed  the  static  torus’ 
residuals  grow  unbounded  to  3,000  km  in  the  same  time.  The  real  value  of  these  results 
cannot  be  known  without  a  direct  comparison  to  a  “full-up”  numerical  propagator  that 
includes  the  same  high-order  gravity  field,  a  sophisticated  atmosphere  model  with  real 
solar  characteristics,  the  complete  solar  and  lunar  ephemerides,  and  solar  radiation 
pressure.  This  recommendation,  left  for  the  next  section,  would  indicate  if  and  how  KAM 
theory  improves  the  current  state  of  the  art.  With  the  present  low-eccentricity  problem,  it 
is  expected  that  the  full-up  integrator  will  be  more  accurate  in  the  short  term  (perhaps  10s 
or  100s  of  meters  in  total  RMS  from  the  true  orbit),  but  the  stochastic  predictions  from 
the  KAM  torus  may  rival  the  propagator  in  the  long  term  as  the  full-up  propagator 
diverges  from  reality. 

Additional  research  will  be  required  to  confirm  these  findings  and  improvements  in 
accuracy  will  be  required  before  KAM  theory  can  be  relied  upon  for  high  priority  assets 
such  as  the  ISS.  This  does  not  preclude  its  immediate  testing  on  debris  orbits  and  low 
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priority  assets  which  could  extend  the  length  of  time  between  reacquisition  by  tracking 
systems. 

5.3  Recommendations  for  Future  Research 

The  first  goal  of  future  efforts  should  be  an  improvement  to  the  automation  process 
for  finding  fundamental  frequencies  from  the  trajectory-following  method.  Once  it  proves 
successful  on  a  handful  of  numerically  integrated  orbits,  a  full  3 -dimensional  matrix  of 
orbits  with  different  eccentricities,  altitudes,  and  inclinations  ought  to  be  assembled  from 
real  satellite  ephemerides  and  batch  processed  to  test  the  reliability  of  the  Fourier 
analysis. 

Another  major  hurdle  that  must  be  overcome  is  the  low  eccentricity  problem 
encountered  here  and  by  Craft.  As  it  relates  to  stochastic  predictions,  it  may  be  possible 
to  improve  the  polynomial  functions  and  reduce  sinusoidal  behavior  in  the  Keplarian  and 
apsidal  coordinates  by  iterating  on  the  polynomial  production  process.  After  the  first 
polynomial  is  generated,  the  Bayes  filtration  process  could  be  repeated  with  initial  torus 
offsets  generated  by  the  previous  polynomial  approximations.  This  procedure  may  be 
iterated  upon  until  the  polynomials  are  smooth. 

The  Bayes  filter  provided  reasonable  results  using  linearized  two-body  dynamics,  but 
it  may  be  worthwhile  to  explore  an  expansion  of  the  equations  of  motion  to  include  the 
moon,  sun,  and  air  drag.  The  current  fitting  process  attempts  to  account  for  all  of  these 
nonconservative  perturbations  by  lumping  them  into  one  set  of  all-inclusive  pseudo¬ 
stochastic  parameters.  Better  results  may  be  obtained  by  isolating  the  parallel  and 
perpendicular  contributions  from  each  source  and  proceeding  to  make  stochastic 
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predictions  from  multiple  sets  of  pseudo-stochastic  parameters.  This  does  not  mean  that 
the  current  dynamics  model  should  be  immediately  discarded,  though.  It  is  quite  possible 
that  some  orbits  (such  as  those  at  higher  LEO  where  air  density  is  reduced)  will  fare  very 
well  with  the  simplified  equations  employed  here. 

As  it  relates  to  the  first  recommendation,  observational  data  ought  to  be  obtained 
from  a  vast  assortment  of  satellite  orbits  and  processed  through  a  Bayesian  filter — similar 
to  the  one  established  here  or  one  improved  by  the  previous  recommendations. 
Predictions  could  then  be  generated  for  each  of  the  satellites  to  determine  the  method’s 
level  of  trustworthiness.  Only  then  will  it  be  certain  that  KAM  theory  can  be  applied 
effortlessly  to  the  multitude  of  objects  that  must  be  tracked  in  earth  orbits.  The  results 
from  this  study  are  simply  not  enough  to  guarantee  the  applicability  of  KAM  theory  for 
operational  use.  This  study  may  also  shed  additional  light  on  the  impacts  of  solar  events. 
Even  though  the  present  study  showed  deteriorations  in  predictions  following  a  major 
solar  storm,  it  may  be  possible  to  compensate  for  deleterious  density  changes  by 
adjusting  the  pseudo-stochastic  parameters  in  a  predictable  manner. 

Finally,  the  results  of  this  study  indicate  that  the  KAM  torus  can  detect  subtle  changes 
in  the  thermosphere  density.  Since  these  changes  are  not  easily  observed  in  the  osculating 
orbital  elements  or  the  native  coordinates,  KAM  theory  may  offer  an  innovative  approach 
for  studying  trends  in  the  upper  atmosphere  from  historical  spacecraft  ephemerides  across 
a  broader  range  of  altitudes  without  the  presence  of  specialized  measurement  devices. 
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Figure  85.  A  silhouette  of  the  Space  Shuttle  Endeavor  (STS-130)  departing  the  ISS  in 
which  the  backdrop  depicts  the  lower  layers  of  the  earth’s  atmosphere.  Credit:  NASA 
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Appendix  A 

A.l  Identification  of  Power  Spectral  Frequency  Combinations  from  0  to  3  rad/TU 


Figure  86.  PSD  plot  identifying  frequency  combinations  from  [0,  0.2]  rad/TU. 


171 


Power, 


Figure  87.  PSD  plot  identifying  frequency  combinations  from  [0.2,  0.4]  rad/TU. 
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Figure  88.  PSD  plot  identifying  frequency  combinations  from  [0.4,  0.6]  rad/TU. 
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Figure  89.  PSD  plot  identifying  frequency  combinations  from  [0.6,  0.8]  rad/TU. 
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Angular  Frequency,  to  (rad/TU) 


Figure  90.  PSD  plot  identifying  frequency  combinations  from  [0.8,  1.0]  rad/TU. 
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Figure  91.  PSD  plot  identifying  frequency  combinations  from  [1.0,  1.2]  rad/TU. 
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Figure  92.  PSD  plot  identifying  frequency  combinations  from  [1.2,  1 .4]  rad/TU. 
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Figure  93.  PSD  plot  identifying  frequency  combinations  from  [1.4,  1.6]  rad/TU. 
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Figure  94.  PSD  plot  identifying  frequency  combinations  from  [1.6,  1.8]  rad/TU. 
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Figure  95.  PSD  plot  identifying  frequency  combinations  from  [1.8,  2.0]  rad/TU. 
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Figure  96.  PSD  plot  identifying  frequency  combinations  from  [2.0,  2.2]  rad/TU. 
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Figure  97.  PSD  plot  identifying  frequency  combinations  from  [2.2,  2.4]  rad/TU. 
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Figure  98.  PSD  plot  identifying  frequency  combinations  from  [2.4,  2.6]  rad/TU. 
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Figure  99.  PSD  plot  identifying  frequency  combinations  from  [2.6,  2.8]  rad/TU. 
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Figure  100.  PSD  plot  identifying  frequency  combinations  from  [2.8,  3.0]  rad/TU. 
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